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SECTION  1.  INTRODUCTION 


This  report  will  describe  progress  made  on  the  STAS  Contract  with  respect  to 
providing  analytical  bounds  and  approximations  to  the  performance  modeling  of 
electro-optical  target  trackers. 


The  major  goal  of  this,  sutdy  was  to  analyze  those  limits  that  one  can  place 
on  an  imaging  tracker  performance  in  terms  of  fundamental  properties  of  the 
sensor  and  the  target  such  as  vehicular  dynamics  and  uncertainty,  signal  to 
noise  of  the  sensor,  and  resolution.  To  accomplish  these  goals  it  was  neces¬ 
sary  to  develop  a  complete  imaging  tracker  system  model  including  a  model  for 
the  relevant  scene,  sensor,  and  observation,  as  well  as  a  model  for  the  tracker 
processor.  It  was  also  necessary  to  establish  those  measures  and  criteria 
by  which  a  tracker's  performance  can  be  gauged  so  that  one  can  produce  from 
these  analyses  a  systematic  design  procedure  applicable  to  the  design,  develop¬ 
ment,  and  specification  of  tracker  systems.  A  number  of  analytic  tools  that 
were  seen  as  potentially  applicable  to  tracker  analysis  were  evaluated  and  a 
subset  of  these  that  were  thought  to  have  particular  application  were  chosen. 
These  tools  were  then  applied  to  the  models  to  produce  bounds  and  approxima¬ 
tions  to  the  actual  performance  one  could  expect  of  a  tracker  in  a  given 
environment. 


The  tasks  that  were  performed  under  this  contract  include: 


1.  Providing  a  simple  analytical  scene  model,  including 
the  effect  of  target  maneuvering,  target  shape 
variation  with  position,  and  clutter  obscuration,  as 
well  as  sensor  degradations  due  to  e.g.,  sensor  blur, 
noise  and  sampling. 

2.  Developing  a  model  that  is  capable  of  bounding  the  optimal 
performance  of  a  tracker.  The  way  this  was  arrived  at  was 
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to  define  In  formal  mathematical  terms  the  function 
of  an  optimal  tracker  and  then  producing  the  equations 
governing  this  function  to  which  there  is  a  uniquely 
defined  exact  solution.  This  solution  to  the  optimal 
tracker  equation  was  then  approximated  and  bounded 
by  looking  at  properties  of  this  solution  and  relating 
It  to  the  performance  parameters  of  the  tracker. 

3.  Establishing  a  set  of  parameters  which  were  capable  of 
gauging  the  tracker  performance.  As  expected,  these 
tracker  performance  parameters  varied  with  the  types  of 
missions  and  scenarios  that  one  would  require  the  use  of 
a  tracker.  A  number  of  these  performance  parameters 
have  been  defined,  and  these  parameters  have  been 
applied  to  sample  tracker  system  sensors  and  the  use 

of  the  measures  have  been  described. 

4.  Taking  all  of  the  key  results  of  the  preceding  three 
tasks  and  verifying  them  by  direct  computer  simulation. 

This  was  performed  by  taking  the  equation  referred  to 
in  (2)  and  numerically  solving  them  and  then  comparing 
the  resultant  solutions  to  our  derived  set  of 
approximations  and  bounds. 

5.  Exploring  the  viability  of  these  techniques  in  actual 
tracker  applications.  There  is  included  in  this  report 
an  example  of  how  these  techniques  can  be  used  for 
approximating,  for  example,  aimpoint  jitter  for  a 
given  precision  tracking  problem. 

Figure  1-1  shows  the  block  diagram  that  will  be  employed  to  model  tracker  systems 
from  a  general  standpoint.  It  starts  off  with  the  scene,  which  can  contain  a 
number  of  targets  and  clutter  objects  as  well  as  various  other  man-made  non-clutter 
objects.  The  scene  Is  then  modeled  as  having  targets  having  a  given  set  of 
dynamics  and  a  given  thermal  signature.  The  effect  of  the  thermal  signature, 
the  projection  of  the  scene,  and  the  emitted  radiation  from  the  relevant  objects 
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In  the  scene  forms  with  the  sensor  model  of  such  effects  as  blur  and  noise, 
the  sensed  Image  produced  by  the  sensor.  This  Image  Is  then  fed  Into  a 
tracker  processor.  The  processor  Is  modeled  as  Implicitly  providing  a  state 
probability  density.  This  may  not  be  an  actual  output  of  the  tracker,  but 
It  Is  treated  In  our  model  as  If  this  density  were  operated  on  to  provide  a  tar¬ 
get  state  estimate  known  as  X.  This  target  state  estimate  can  then  be  used  to 
drive  the  servomechanism  which  can  control  a  set  of  gimbals,  which  can  then  be 
fed  back  through  a  gimbal  pickoff  to  feed  a  set  of  derived  angles  into  the 
tracker.  There  are  two  types  of  trackers  one  can  describe,  known  as  either 
open  loop  or  closed  loop.  In  the  closed  loop  case,  the  servomechanism  directly 
drives  the  gimbals  which  controls  the  pointing  of  the  sensor  and  which  can 
feed  information  to  the  tracker  as  to  the  new  sensor  line  of  sight  position. 

The  open  loop  case  is  tracking  "in  raster"  and  basically  treating  the  sensor 
as  if  it  were  stationary.  It  can  be  seen  that  the  transformation  of  an  open 
loop  to  a  closed  loop  tracker  Invloves  Including  a  certain  amount  of  pickoff 
noise  and  analyzing  the  relevant  servomechanisms.  Although  these  problems 
are  non-trivlal,  they  are  beyond  the  scope  of  this  particular  report  since  it 
is  felt  that  the  analysis  of  servomechanisms  is  treated  in  sufficient  detail 
in  many  other  publications  and  reports.  The  open  loop  case  will  therefore  be 
treated  in  this  report  although  the  generalization  of  our  formulation  to  the 
closed  loop  case  will  be  apparent. 

To  discuss  In  more  detail  the  modeling  of  a  tracker,  refer  to  Figure  1-2, 
where  a  series  of  sensed  images  is  represented  as  the  vector  function 
■y(t)  and  the  tracker,  as  previously  discussed,  produces  a  probability 
estimate  of  the  statelet),  P0T(t),t).  T(t)  reoresents  the  state  of  the  tar¬ 
gets  and  of  the  scene  that  one  Is  tracking.  The  state  includes  such  things  as 
target  position,  velocity,  and  acceleration,  in  addition  to  clutter  positions 
and  obscuration  co-efficients.  In  general,  the  state  vector  changes  as  a 
function  of  time.  The  resultant  state  probability  density  can  then  be  operated  on 
by  an  arbitrary  functional  operator  producing  a  state  estimate.  Examples  of 
the  types  of  estimators  one  can  use  are,  e.g.,  a  minimum  mean  square  estimate  which 
would  take  the  mean  of  the  underlying  probability  distributions,  or  a  maximum 
a  priori  estimator  which  could  pick  the  peak  of  the  probability  density.  One 
could  also  define  a  minimum  error  probability  estimator  or  a  mini-max  type 
function  as  well  In  crder  to  produce  an  estimate  that  minimizes  the  criteria 
of  parties  ar  relc  .e  In  the  scenario  of  interest. 
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FIGURE  1-2.  TRACKER  MODELING 


This  report  will  be  broken  into  eight  sections.  After  the  introduction, 
the  second  section  describes  the  analytic  techniques  that  have 
potential  value  for  the  modeling  of  tracker  systems  in  various  environments. 
These  techniques  include  the  ways  of  analyzing  the  stochastic  differential 
equations  that  are  central  to  the  definition  of  the  tracker  problem.  The 
third  section  will  describe  the  target  model  including  target  dynamics, 
maneuvering,  and  shape.  The  fourth  section  is  related  to  scene  modeling 
and  will  describe  the  techniques  for  taking  the  scene  state,  (which  includes 
the  target  state  as  well  as  the  clutter  and  obscuration  states)  and  producing 
from  the  target,  clutter  and  scensor  models,  an  expression  for  what  the 
sensed  IR  Scene  will  be  in  terms  of  the  sampled  output.  The  fifth  section 
deals  with  bounds  and  approximations  and  will  describe  various  techniques  for 
bounding  and  approximating  the  performance  of  the  optimal  filter  operating 
on  the  video  provided  by  the  scene  model.  Following  this,  the  sixth  section, 
titled  Simulation,  will  describe  the  computer  program  that  has  been  developed 
to  directly  simulate  the  exact  solution  to  arbitrary  non-linear  stochastic 
differential  equation,  how  it  has  been  applied  to  the  problems  of  interest, 
and  some  of  the  results  that  have  been  obtained.  In  the  seventh  section, 
system  applications  of  the  techniques  described  in  this  report  will  be 
discussed.  Arid  finally,  a  summary  and  conclusion  of  the  contents  of  the  report 
will  be  given. 
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.  SECTION  2.  BASIC  ANALYTIC  TECHNIQUES 
In  this  section,  the  basic  analytic  techniques  that  can  be  used  to  model, 
bound  and  approximate  the  performance  of  FLIR  trackers  are  summarized.  First, 
the  fundamental  nature  of  the  tracker  modeling  problem  is  discussed,  along 
with  a  definition  of  the  interpretation  of  the  meaning  of  the  word  'tracker' 
used  in  this  report.  A  set  of  representative  tracker  evaluation  parameters 
and  normalizing  system  parameters  are  then  given.  This  leads  to  a  formulation 
of  the  problem  in  terms  of  stochastic  differential  equations.  The  basic 
characteristics  of  these  types  of  equations  are  summarized  along  with  the  two 
types  of  disturbance  processes,  Wiener  and  Poisson,  applicable  to  the  tracker 
problem.  Then  the  extent  to  which  a  solution  can  be  found  is  discussed  in 
terms  of  Kushner's  equation  and  its  reduction  to  the  Kalman  filter  in  the 
linear  case.  Finally,  approaches  for  computing  performance  bounds  and  approxi¬ 
mations  for  certain  classes  of  optimal  nonlinear  trackers  are  presented 
and  discussed. 


2-1 


2.1  TRACKER  MODELING  PROBLEM 

Before  the  fundamental  nature  of  the  tracker  modeling  problem  can  be 
defined,  it  is  necessary  to  define  the  meaning  of  the  word  'tracker',  in  the  context 
of  this  analysis  study.  In  general  terms,  trackers  can  be  classified  as  either 
open-loop  or  closed-loop.  With  closed  loop  trackers,  the  sensor  is  driven  to 
maintain  the  target  in  a  fixed  position  in  the  sensor  field  of  view.  In  an 
open-loop  tracker,  the  location  of  the  target  in  the  FOV  is  tracked  as  it  moves 
through  the  FOV.  In  this  report,  the  results  will  be  derived  for  an  open-loop 
type  tracker  operating  on  FLIR  video  of  a  fixed,  noisy,  cluttered  scene  with  a 
target  moving  through  the  scene.  This  restriction  is  made  for  simplicity  of 
development,  with  no  real  loss  of  generality  since  the  results  and  techniques 
developed  can  easily  be  extended. 

The  basic  tracker  modeling  problem  is  to  develop  a  mathematical  model  of 
tracker  performance  that  allows  extraction  of  figures  of  merit  for  use  in  comparing 
the  relative  performance  of  various  tracker  implementations  with  bounds  on  optimal 
tracker  performance  for  various  applications.  The  model  must  be  broad  enough  to  " 
encompass  the  analysis  of  relatively  simple  centroid  trackers,  correlation" 
trackers,  and  advanced,  intelligent  target  trackers  now  under  development  (cf  ref  (27)) 
The  latter  class  of  trackers  use  modern  image  processing  and  pattern  recognition 
techniques  to  initially  acquire  targets  and  to  predict  or  mitigate  occlusion 
effects  of  clutter  also  present  in  the  scene.  In  this  sense,  this  problem  is  a 
natural  extension  of  a  previous  study  on  autocuer  modeling  analysis  (ref  (30)). 

The  extension  is  basically  to  a  time  history  of  observations  being  used  to  predict 
future  time  history. 
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In  order  to  illustrate  the  scope  of  application  and  its  affect  on  appro¬ 
priate  figures  of  merit,  two  examples  of  the  use  of  FUR  trackers  are  pertinent: 
High  Energy  Lasers  (HEL),  and  advanced  tank  fire  control. 

In  the  HEL  application,  the  requirement  is  to  provide  automatic  aimpoint 
selection  and  maintenance  for  the  HEL  beam.  Performing  this  function  requires,  in 
turn,  highly  accurate  position  tracking  to  maintain  the  beam  on  one  spot;  and 
target  shape  and  orientation  tracking  to  find  the  proper  vulnerable  aimpoint. 
Accurate  target  dynamics  prediction  is  important  to  the  extent  that  it 
affects  the  tracker  stability. 

In  the  advanced  tank  fire  control  application,  the  requirements  will  be 
to  automatically  find,  identify,  and  track  targets;  while  providing  capability 
to  hit  maneuvering  targets  and  automatically  position  the  gun  to  hit  the  target. 
Here  the  critical  requirement  is  to  accurately  predict  the  target  location  in 
the  future  at  the  intercept  point,  typically  one  to  two  seconds.  Meeting  this 
requirement  requires  accurate  modeling  and  estimation  of  target  dynamics.  The 
effect  of  not  properly  modeling  target  dynamics  is  Illustrated  in  (32), 
where  it  is  shown  that  an  extended  Kalman  filter  implementation  of  a  tracker  could 
not  follow  an  abruptly  maneuvering  target  in  the  air-air  scenario. 

It  is  apparent  from  these  two  examples  that  different  figures  of  merit 
apply  to  different  applications.  The  appropriate  list  of  figures  of  merit, 
however,  can  be  easily  identified: 

1)  Expected  RMS  tracker  error 

2)  Stability  (probability  of  losing  lock) 

3)  Prediction  accuracy 

4)  Robustness  (model  sensitivity) 

These  figures  of  merit  will  be  functions  of  the  conditions  under  which 
the  tracker  will  operate: 
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1)  Signal  to  noise  ratio 

2)  Clutter  density 

3)  Target  maneuver  parameters 

4)  Target  size 

5)  Degree  of  non-linearity 

In  order  to  relate  the  figures  of  merit  to  these  parameters,  a  model  of 
the  tracker  performance  must  be  developed.  This  can  be  done  formally  as  follows: 
In  section  3,  a  3-D  target  dynamics  model  is  developed  which  has  the 

form 

dx“  =  f(T,  t)dt  +  diT  (2-1) 


where  x  is  the  target  dynamics  state  vector  (augmented  by  clutter  states 
as  described  in  section  4)  and  w  is  an  appropriate  random  process. 

Also  in  section  4,  a  2-dimensional  projection  of  a  target  at  state  x  on 
the  sensor  image  plane  is  shown  to  have  the  general  form 


9(u,  v)  =  hQ(u,  v,  x) 
where  u,  v  are  the  sensor  coordinates 


(2-2) 


In  section  4,  a  measurement  model  is  shown  to  obey  an  equation  of  the  form 


I(u,  v) 


(2-3) 


where  I(u,v)  is  the  intensity  function  at  (u,v), 
th 

y^  is  the  i  component  of  vector  y  and  y  is  defined  by 

dy-  =  h(x)dt  +  dz  (2-4) 


h(x) 
with  u^ ,  v^ 


=  [h0(ur  v]t  x),  ho(u2,  v2,  x),  ho(u3,  v3,  x)  ...  ]T 
being  the  u,v  coordinates  of  the  ith  pixel. 
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Equations  (2-1)  and  (2-4)  formally  define  the  system  in  terms  of  stochastic 
differential  equations,  where  the  tracker  synthesis  problem  is  to  find  a  near 
optimal  estimate  of  x  based  on  measurements  y  and  the  performance  estimation  problem 
is  to  estimate  or  bound  the  performance  of  such  estimates. 

2.2  STOCHASTIC  DIFFERENTIAL  EQUATIONS 

Analysis  of  the  tracking  problem  in  terms  of  stochastic  differential  equa¬ 
tions  requires  the  introduction  of  special  techniques,  since  the  differential  of 
the  random  processes  involved  does  not  exist  in  the  sense  of  ordinary  differentials. 
The  formulation  that  provides  meaning  to  these  special  differentials  is  generally 
referred  to  as  It6  calculus  and  the  interested  reader  is  referred  to  any  of  the 
general  texts  on  this  subject  listed  in  the  bibliography,  e.g.,  McGarty  (2). 

Here,  the  basic  concepts  useful  in  tracker  analysis  are  summarized  to  provide 
guidance  for  further  study. 

The  first  concept  is  that  of  independent  increments.  Basically,  this  con¬ 
cept  means  that  increments  of  the  process  over  non-overlapping  (time)  intervals 
are  independent.  Formally,  the  definition  is  as  follows: 

Let  for  t .  <  t .  s  t  <  t 
J  i  m  n 

V  =  x(t.)  -  x ( t j )  (2-6) 

"  ■  x(t„)  -  *<*»> 

Then  x ( t )  is  an  independent  increment  process  if  v  and  u  are  independent  random 
variables. 

The  importance  of  this  concept  lies  in  the  fact  that  two  well  known  random 
processes  that  are  directly  applicable  to  modeling  the  tracker  problem  satisfy 
these  conditions.  These  are  the  Wiener  and  Poisson  processes. 

The  Wiener  process  is  useful  both  for  modeling  noise  in  measurements 
model  and  disturbances  in  the  dynamical  model.  This  is  the  most  common  process 
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used  in  estimation  theory,  and  has  the  following  formal  definition: 

x ( t )  is  a  Wiener  process  if  it  has  independent  increments,  and 
x (  0)  =  0 


p[(x(t)  -  x( S ) )  ?  X] 


oy  2 


±—f 

Tr(t-s)  J 


e-n2/(2o2(t-s))dr| 


(2-7) 


for 

t  >  s 

As  a  consequence,  then,  this  process  has  the  following  important  proper¬ 
ties  : 

E[x ( t) ]  =  0 

E[x2(t)]  =  a2t 

E[x(t)  x(s)J  =  min  (t,  s) 


and,  if 

AX  =  X ( t  +  At)  -  x(t)  (2-8) 

then 

EfAx2]  =  crZAt  (2-9) 

The  Poisson  process  (2)  '  is  useful  in  modeling  step  changes  in  the 

system  state  caused  by  execution  of  maneuvers  at  random  times.  The  formal 

definition  of  a  Poisson  process  is  as  follows: 

x(t)  is  a  Poisson  process  if  it  has  independent  increments,  the  sample 

* 

function  x(t,w)  is  a  step  function  increasing  with  jump  1,  and: 

*  A  Poisson  process  can  be  generalized  to  have  several  jump  values  a<  each  with 
a  characteristic  probability  of  occurrence  Pg 
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x(0)  =  0 


p [ x ( t )  -  x(s)  2  k] 


-x| t-s I  (|t-sl)kxk 
e  k! 


where  >  >  0 


(2-10) 


Also 


x(t-)  t  x(t)  =  x(t+) 
i.e.  continuous  on  the  right. 


(2-11) 


The  special  properties  of  these  processes  lead  to  a  requirement  to 
extend  the  definitions  of  differential  and  integral  Stochastic  processes. 

As  shown  above,  e.g.,  the  second  moment  of  the  Wiener  process  is  of  the  order  dt, 

2 

not  dt  so  that  higher  order  terms  must  be  retained  in  any  expansions.  This 
property  also  results  in  the  process  not  being  of  bounded  variation,  so  differ¬ 
entials  do  not  exist  in  the  ordinary  sense.  The  solution  to  this  problem  was 
developed  by  Ito  and  is  illustrated  here  by  the  Ito  rule  for  differentials,  pre¬ 
sented  here  in  the  scalar  form. 

Let 


where 


and 


y  =  f(x)  (2-12) 

dx  =  g(x)  dt  +  dw  (2-13) 

E[dw  dw]  =  Wdt  (2-14) 
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The  problem  is  to  find  the  differential  dy.  The  approach  is  to  perform  an 
expansion  through  second  order  terms. 

Then 

AY  =  f ( x  +  ax)  -  f(x) 


-  df  .  .  1  / 

'  di  4x  *  2  ( 


Ax 


(2-15) 


5  +^t?(9(X)At  +AW)‘ 


Retaining  only  first  order  terms  in  At  results  in 


+  j  ~  W  At 
c  dx^ 


(2-16) 


which,  by  taking  limits,  results  in 


■  (&)“’- 4(0)“dt 


(2-17) 


Note  that  dw  was  replaced  by  its  expectation,  since  with  probability  1, 

„t 


/ 


(dw)^  Wt 


for  any  finite  t. 

With  this  understanding  of  stochastic  differentials,  the  stochastic 
differential  equation  formulation  of  the  tracker  problem  given  formally  in  sec¬ 
tion  2.1  has  meaning.  The  problem,  then,  is  from  observing  the  measurement  y 
estimate  the  probability  density  of  the  state  vector  x.  Then  the  minimum  mean 
square  error  (MMSE),  or  maximum  aposteriori  (MAP)  estimate  of  x  can  be  deter¬ 
mined.  In  addition,  from  the  probability  density  of  x,  the  error  can  be 
estimated. 
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In  developing  a  model  of  a  system  in  terms  of  stochastic  differential 
equations,  an  additional  effect  needs  to  be  considered.  Suppose  the  system 
can  be  represented  in  the  form 


dx  =  f(x,t)  dt  +  g(x,t)  dnB 


(2-18) 


where  nD(t)  is  integrated  white  noise  of  bandwidth  B  and  w(t)  is  a  Wiener 

D 

process  satisfying 


/ 


limit  /  (dn0  -  dw) 

B-k»  '  a 


(2-19) 


Then 


/■ 


(dx  -  f(x,t)  dt)  = 


limit  /  (g(x,t)  |f  dx)  dng 
imit  /  (g(x,t)  +  j  g(x,t)  dng)  dn 

B-*»  J 

imit  /  (g(x,t)  +  |  |f  g(x,t)  dw)  dw 

B-h»  J  c 

-  J g(dw  -  dng)  -  \  j ff  9  (dw2  -  dng2) 
=  j (9(x,t)  +  \  | J  gdw)  dw  -  j  J ||  gdt 


=  1 


=  li 
B-x» 


(2-20) 


The  last  term  is  referred  to  as  the  Wong-Zakai  correction,  so  in  the  limit 
as  b  -*■  *  ,  x  obeys  the  differential  equation 


dx  =  [f (x.t)  -  ^  9  ff]  dt  +  9(x,t) 


(2-21 


***  ' 


2.3  SOLUTION  OF  THE  ESTIMATION  PROBLEM 

As  stated  in  section  2.2,  the  basic  estimation  problem  is  to  find  the 
probability  density  function  P  of  the  state"?  given  the  observation  or 
measurement  ' y .  There  are  two  useful  approaches.  Kushner  has  developed  a 
complete  solution  for  P  in  terms  of  a  stochastic  partial  differential  equa¬ 
tion  (ref  (43))  and  the  Bucy  representation  theorem  can  be  used  for 
discrete  simulation  of  results  on  a  digital  computer.  The  Kushner  solution 
is  given  here  for  an  estimation  problem  with  generalized  Poisson  process  in 
addition  to  Wiener  processes.  The  formal  definition  of  the  basis  estimation 
problem  is  given  in  Table  2-1,  along  with  Kushner's  solution  for  the 
probability  density  function  P. 

This  equation  is  a  nonlinear  stochastic  partial  differential  equation 
and  has  no  general  solution.  Thus  the  problem  is  in  finding  tractable  analytic 
means  for  dealing  with  it.  A  simple  graphical  illustration  of  the  terms  in  the 
Kushner  equation,  as  shown  in  Table  2-2,  helps  to  provide  an  insight  into  the 
nature  of  this  problem.  The  first  term,  A  (y),  is  a  matched  filter  term 


TABLE  2-1.  FORMULATION  AND  FORMAL  SOLUTION  OF  THE 
OPTIMAL  ESTIMATION  PROBLEM  USING 
STOCHASTIC  DIFFERENTIAL  EQUATIONS 


Let  dx  *  f(x,t)dt  +  dw  +  dn 
dy  ■  h(x,t)dt  +  dz 


w, 

dn 


z  Wiener  process  with  E[dw  dw]  *  Wdt,  E[dz  dz]  «  Zdt 
Poisson  process  with  rate  X  and  transition  density  Pa 


i 


Let  P  =  P(x,t)  represent  probability  density  of  x(t)  at  time  t 


Then 


dP  *  A(y)  +  B(P)dt  +  C(P)dt 
Where 


A(y)  -  P(dy  -  Hdt)Z-1  (h  -  H)  "Matched  Filter" 

B(P>  '  Esf;  <fip>  +  <w>  "Orlft-DiffusW 


C(P)  -  [-P  +  P  @  P]  "Abrupt  Change" 

1 


TABLE  2-2.  GRAPHICAL  ILLUSTRATION  OF  TERMS  IN  THE  KUSHNER  EQUATION 


modified  by  the  probability  function  P(x,t).  This  is  illustrated  in  Table  2-2 
for  two  different  measurements  dy  and  dy^.  The  roll-off,  as  the  observation 
deviates  from  H,  the  expected  value  of  h  (x),  is  due  to  the  P  (x,t) 
modification. 

The  second  term  B(p)  is  a  combination  drift  and  diffusion.  This  is  a 
noise  effect  and  is  represented  by  a  second  order  partial  derivative,  since 
it  is  a  relative  effect.  For  example,  as  shown  in  the  table,  the  net  change 
in  the  region  is  the  resultant  of  the  diffusion  from  A,  into  A^  minus  the 
diffusion  from  A^  into  A^  and  similarly  for  regions  ^  and  A^.  The  drift 
represents  the  fact  that  f  (x,t)  will  cause  the  probability  distribution  to 
have  a  net  shift  as  shown  in  the  figure. 

The  third  term  is  the  Poisson  jump  terms.  Referring,  again,  to  the 
table,  the  arrows  represent  the  probability  of  a  jump  of  magnitude  equal  to  the 
x  coordinate  of  the  arrow.  This  is  convolved  with  the  old  probability  density 
(solid  line),  resulting  in  a  new  probability  density  (one  of  the  dashed  lines). 
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This  is  then  modified  by  the  probability  \±  that  a  jump  will  occur  in  that  time 
interval ,  resulting  in  the  Poisson  jump  contribution  to  the  Kushner  equation* 

In  addition  to  the  continuous  time  solution  represented  by  the  Kushner 
partial  differential  equation,  there  is  a  discrete  time,  integral  equation 
formulation  (see  Table  2-3) ,  and  its  solution  is  equivalent  to 
the  solution  of  Kushner's  equation.  This  discrete  solution  is  more  useful  for 
digital  simulation,  which  by  it's  very  nature  is  discrete.  This  solution  is 
based  on  the  representation  theorem. 

Representation  Theorem 

The  representation  theorem  first  proven  by  Bucy,  expresses  the  conditional 
densityin  an  integral  which  could  be  used  to  derive  results  which  would  be 
difficult  or  impossible  using  Kushner’s  equations.  In  fact,  Kushner’s 
equations  can  be  derived  from  this  theorem.  Using  the  representation 
theorem,  we  get  for  Tt  and  defined  as  in  section  2,  and  for 
X  =  0 

Px(u.t|y(T).  T.t)  =  t)  tf-22) 

X  E[eHt]  x 

with  EC  3  denoting  expectation  and 

Ht  =  J  tiT(x(n),n)  Z_1(dy(n)  -  ^  h(x(n),  n)  dn)  (2-23) 

-oo 

The  compactness  and  elegance  of  equations  (2-22)  and  (2-23)  allows  this 
expression  for  the  state  probability  density  to  be  used  effectively  in 
proofs  related  to  the  stochastic  filtering  problem  (see  (2)). 


t 
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x(K)  =  f(x(K-1))+  Wi(K  -  1) 

z(K)  =  h(x(K  -  1))  +  W,(K  -  1)  W1‘  W2  GAUSSIAN 


o 


LL 

O 

□ 

O 

O 
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I! 


\o  o 

-I  -J 


R  =  E[W2(K)  W2(K)1 

Q  =  E[W-|(K)  W-j(K)] 

N(A,  B)  =  EXP  [(-1/2)  •  (ATB‘1A)] 


Although  neither  the  integral  nor  differential  form  of  the  solution  can  be 
solved  in  general,  there  are  some  special  cases  where  a  solution  can  be 
obtained.  For  the  case  of  linear  dynamics  and  linear  measurements  with 
Gaussian  noise,  it  reduces  to  the  familiar  Kalman-Bucy  filter.  Similarly, 
many  non-linear  systems  may  be  approximated  by  a  linearization  of  the  form 

?(*,t)  =  *ft**,t)  +  A(t)0?-**) 

A  *  (2-24) 

tor.t)  =T(t*,t)  +  c(t  )(*-**) 

and  the  equation  will  reduce  to  the  extended  Kalman-Bucy  filter  for  estimating 
Explicit  formulations  for  these  cases  can  be  found  in  (  3  ). 
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2.4  PERFORMANCE  BOUNDS  AND  APPROXIMATIONS 

For  non-linear  problems,  a  closed  form  solution  of  the  Kushner  equation 
cannot  be  produced,  in  general.  The  problems  of  calculation  of  the  optimal 
non-linear  estimator,  or  even  the  calculation  of  the  performance  of  simple  sub- 
optimal  estimators,  are  computationally  intractable,  since  both  would  require 
infinite  dimensional  calculation  (223-  The  principal  problem  in  calculating 
performance  is  that  the  propagation  of  the  second  moment  includes  not  only 
the  past  time  history  of  the  second  moment,  but  also  of  all  higher  moments. 

This  results  in  an  infinite  chain  of  higher  order  moments  that  must  be  propa¬ 
gated  in  order  to  establish  mean  square  error  performance.  This  behavior  is 
in  dramatic  contrast  to  the  situation  with  linear  systems  where  the  second 
moment,  or  mean  square,  error  propagates  by  a  Riccati  type  differential 
equation,  and  higher  order  moments  have  no  effect  on  performance  evaluation. 

In  order  to  circumvent  the  problems  inherent  in  the  non-linear  tracker  esti¬ 
mation  problem,  two  related  analytic  techniques  that  are  both  practical  and 
mathematically  rigorous  can  be  used:  performance  bounds  and  approximation  by 
linearization. 

Performance  bounds  are  generated  by  making  controlled  approximations 
that  guarantee  that  the  resulting  performance  estimate  is  either  less  than 
(lower  bound)  or  greater  than  (upper  bound)  the  performance  of  the  optimal 
estimator,  or  for  any  candidate  suboptimal  estimates.  Such  bounds  can 
always  be  found,  but  the  analytical  trick  is  to  make  the  approximations 
judiciously  so  that  the  resulting  bounds  are  in  some  sense  "tight"  and  non¬ 
trivial.  A  tight  bound  can  be  assured  if  the  difference  between  the  upper 
and  lower  bounds  is  small. 
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Good  performance  bounds  can  be  extremely  valuable  to  the  system  designer. 
If  the  lower  bound  of  a  candidate  suboptimal  design  is  close  to  the  lower 
bound  on  optimal  performance,  he  knows  he  has  done  his  job  well  and  further 
refinements  or  improvements  may  not  be  cost  effective  or  may  be  futile. 
Alternately,  if  the  required  performance  is  lower  than  the  lower  bound  on  the 
optimal  system,  the  futility  of  the  system  design  problem  becomes  apparent. 
Upper  bounds  are  useful  in  that  if  the  upper  bound  of  a  candidate  design  is 
lower  then  the  system  requirements,  then  it  is  not  necessary  to  seek  further 
improvements  in  design.  The  performance  bounds  also  provide  the  all-important 
functional  relationship  between  system  parameters  and  performance.  This  allows 
the  designer  to  perform  tradeoffs  on  the  parameter  requirements  to  insure  a 
cost  effective  system.  For  the  specific  application  to  tracker  systems,  these 
bounds  can  be  obtained  for  both  the  mean  square  tracking  accuracy  and  the 
probability  of  losing  lock,  both  important  figures  of  merit  for  trackers. 

The  second  useful  analytic  technique  is  to  approximate  the  non-linear 
estimation  problem  with  a  linear  problem,  through  a  linearization  technique. 
This  approach  also  provides  a  functional  relationship  between  system  para¬ 
meters  and  performance,  and  provides  an  estimate  of  the  actual  performance. 

This  technique  can  be  used  in  conjunction  with  the  performance  bound  techni¬ 
que.  If  the  bounds  are  tight  and  the  performance  of  the  extended  Kalman-Bucy 
filter  for  the  linearized  system  falls  between  the  upper  and  lower  bounds, 
then  the  performance  of  the  linearized  system  must  be  close  to  the  optimal 
non-linear  estimator.  This  is  often  the  case  and  then  the  linearization 
approach  provides  the  design  solution.  Linearization  is  also  important  in 
that  it  is  often  used  as  a  step  in  computing  performance  bounds,  as  shown 
below.  An  example  of  the  use  of  linearization  is  given  in  section  5. 
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Two  general  techniques  have  been  identified  and  investigated  for  com¬ 
puting  mean  square  error  performance  bounds  for  the  non-linear  estimation 
problems.  Each  of  these  approaches  relates  performance  bounds  from  the 
optimal  non-linear  system  to  performance  of  related  linear  systems. 

The  first  approach,  the  Bobrovsky-Zakai  bound  (23)*  is  based  on  an 
infinite-dimensional  extension  of  the  Van  Trees  version  of  the  Cramer-Rao 
bound  (42).  This  extended  bound  is  used  to  define  a  linear  system  related 
to  the  original  non-linear  system,  but  not  a  direct  linearization,  that  has 
the  property  that  the  optimal  solution  of  the  linear  system  has  an  MMSE 
that  forms  a  lower  bound  for  the  optimal  non-linear  estimation  error.  An 
example  of  this  approach  for  a  simple  one-dimensional  problem  is  presented 
in  section  5,  with  details  of  the  derivation  in  appendix  E. 

The  second  approach  is  termed  the  cone-bounded  approach,  and  it  provides 
both  upper  and  lower  bounds  on  the  performance  of  the  optimal  non-linear  esti¬ 
mator  in  terms  of  the  performance  of  the  associated  linearized  system  and 
certain  non-linearity  parameters  that  are  representative  of  the  degree  of  the 
non-linearities.  These  bounds  apply  to  a  restricted  class  of  non-linear  sys¬ 
tems  that  satisfy  a  uniform  Lipschitz  condition  of  the  form 


|  | f (x+Ax,t)  -  f(x,t)  -  A(t)  Ax | |  £  a(t)|j  Ax | | 
j |h(jd-Ax,t)  -  h(x,t)  -  C(t)  Axj  |  ^  e(t) | |  Ax| j 


(2-28) 


Although  this  approach  is  restricted  to  systems  that  satisfy  (2-28),  the 
approach  is  still  useful  since  similar  Lipschitz  conditions  are  usually  needed 
to  assure  unique  solutions  to  ordinary  non-linear  differential  equations. 

Since  physical  systems  are  generally  of  this  type,  this  approach  will  be  useful 
in  practical  applications. 
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The  result  for  the  lower  bound  has  been  extended  in  this  project  to  give 


the  following  result.  Let  the  linearized  system  be 


dx  =  Axdt  +  dw 
dy  =  Cxdt  +  dZ 


(2-29) 


and  let  P  be  the  error  covariance  matrix  associated  with  the  optimal  filter  for 
this  system.  Then,  the  lower  bound  for  the  non-linear  system  covariance,  Q, 
is  given  by 


Q  >  (1  -  r(t) )  P(t)  (2-30) 

where  r(t)  is  the  unique  non-negative  root  of  the  equation 

(1  -  r(t))  er(t)  -  e-d(t)  (2-31) 

and 

d(t)  -  f  [a2(S)  W_1(S)+  c2(S)  Z-1(S)]  T.[p(S)]  uS  (2-32) 

U 

A  A 

where  T  (P)  ■=  Trace  of  P. 
r 

This  is  a  recent  improvement  ^44 j  over  the  original  result  ^22j  where  d(t)  was 
computed  using  the  trace  of  the  upper  bound  on  performance.  Thus,  the  lower 
bound  is  both  tighter  and  much  simpler  to  evaluate. 

Note  chat  if  a  and  6  are  defined  as 


-t  2  -1 

-  *n[a(t)]  -  J  a  (S)  W  A(S)  Tr  P(S)  dS 

(2-33) 

-  fcn[6(t)]  -  f  c2(S)  Z-1(S)  Tr  P(S)  dS 


then  r(t)  satisfies 

(l-r(t))  er<t)  -  a(t)  6(t)  (2-34) 

Thus,  a  and  8  are  measurements  of  the  degree  of  non-linearity  of  the  dynamics 
and  measurement  equations,  respectively. 
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The  principal  problem  in  using  this  approach  is  in  maintaining  a  tight 
bound  in  the  Lipschitz  condition  (equation  2-28).  If  the  a(t)  or  c(t)  func¬ 
tions  are  too  large,  then  both  the  upper  and  lower  bounds  will  quickly 
degenerate  into  trivial  bounds.  This  is  a  consequence  of  equation  (2-32), 
which  shows  that  if  these  functions  are  too  large,  then  d(t)  will  become  cone 
bounded,  forcing  r(t)  to  approach  1.  The  lower  bound  then,  e.g.,  becomes  0 
from  equation  (2-30). 
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2.5  BOUNDS  ON  PROBABILITY  OF  LOSING  TRACK 

The  techniques  discussed  in  Section  2.4  applied  to  finding  bounds  on  the 
accuracy  of  tracking,  and  as  such,  are  not  applicable  to  the  problem  where  the 
tracker  has  completely  lost  track.  Thus,  an  additional  performance  parameter 
is  needed.  This  parameter  is  the  probability  of  losing  track  in  some  time 
interval.  This  is  analogous  to  the  situation  in  the  ATAC  Autocuer  Modeling 
Analysis  OOJ  where,  in  addition  to  accuracy  in  locating  an  object,  the  proba¬ 
bility  of  totally  missing  the  object,  called  the  probability  of  anomaly,  was 
needed  to  provide  complete  performance  analysis. 

The  probability  of  losing  track  is  relative  to  the  concept  of  non-linear 
stability.  This  concept  is  illustrated  in  Figure  2-1.  Basically,  the  concept 
involves  answering  the  question:  Given  a  non-linear  system  with  random  inputs 
and  outputs,  what  is  the  bound  on  the  probability  that  the  random  process  output 
is  greater  than  some  threshold  in  some  interval,  given  its  initial  value?  Here, 
an  approach  to  finding  an  upper  bound  to  the  probability  of  losing  track  is 
presented. 

A  more  in-depth  discussion  of  this  method  can  be  found  in  Kushner  £0 . 

The  basic  method  can  be  applied  to  any  stochastic  process  which  satisfies  an 
equation  of  the  form 


dX  =  ?(X,t)dt  +  o(X,t)dw 
where  X  is  a  n-dimensional  state  vector, 
o(Z,t)  is  a  nxm  matrix 

w,  a  m-dimensional  normalized  Wiener  process. 


(2-35) 


r 


For  applications  to  the  tracker  problem,  the  state  variable  in  (2-35)  will 
generally  be  taken  to  be  the  error  process 


e  =  x  -  x 

where  x  is  the  target  location  and  x  is  the  target  location  estimate.  In  gener¬ 
al,  it  will  not  be  possible  to  write  the  error  process  in  the  state  equation  form 
de  =  "f(e,t)dt  +  o(e,t)dw  .  (2-36) 

However,  a  linearization  of 
de  =  dx  -  dx 

about  the  state  estimate  x  will  result  in  an  equation  of  form  (2-36).  This 
is  necessary  in  order  to  apply  the  methods  in  £lj. 

For  the  tracking  problem,  one  will  be  interested  in  the  probability 


that 


This  probability  is  bounded  via  methods  in  [  l]  by  first  applying  the  differential 


generator  A  to 


where  A  -  Z  ^(e.t)  +  ££(cr  aT) 


ij  ^ 


f  ^  is  the  itl1  element  of  the  vector  f  and  (a  a^) ^  is  the  i,j  element 

of  the  matrix  a  oT.  The  differential  generator  A  can  be  treated  as  the 

expected  value  of  the  time  derivations  of  its  argument.  If  positive  y  and 

*  can  be  found  such  that 
t 

A( j e | |<  -p  | |e| |  +  $(t) , 
then  the  upper  bound  is  given  by 


supj|e(s)!|>m  |  J"e( o) )  )=k)j<  -  M  exp  f-  I  f 

o<5^t  /  \m  l 
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In  effect  the  y  can  be  looked  at  as  a  degenerative  feedback  term  on 
the  error  norm  ||e|j  and<j)(t)can  be  modeled  as  the  degree  to  which  the  tracker 
departs  from  an  Idealized  one  pole  feedback  system. 

The  overall  approach  to  using  this  concept  of  non-linear  stability  to 
compute  the  probability  of  losing  lock  is  illustrated  in  Figure  2-2.  The 
basic  problem  is  to  find  a  bound  on  the  random  process  e,  which  is  the  differ¬ 
ence  between  the  estimate  of  the  state  and  the  true  state.  Having  a  bound  on 
this  error  provides  a  means  for  bounding  the  probability  of  losing  lock,  and 
the  method  discussed  above  provides  one  means  of  doing  this. 

The  key  is  an  upper  bound  on  the  probability  that  in  the  time  interval 
(0,t)  the  error  will  exceed  a  specified  value,  given  that  the  error  is  known 
at  time  0.  The  procedures  applies  generally  to  any  suitable  function  G(e),  not 
just  the  norm.  Selecting  this  function  (similar  to  selecting  a  Liaponov  function 
for  other  stability  applications)  is  a  key  step  in  the  procedure  as  is  selecting 
<J>  which  can  be  quite  arbitrary.  For  linear  systems,  some  insight  into  the 
nature  of  the  problem  must  be  used  since  these  parameters  basically  relate  to 
the  degree  of  non-linearity.  The  procedure  can  always  be  applied,  some  y  and  <J>(t) 
can  always  be  found,  however,  if  they  are  not  selected  carefully,  the  bound 

may  prove  to  be  trivial.  Note  that  since  the  process  is  applied  to  the  error 
in  the  estimation,  all  normalizing  factors  previously  identified  in  the  per¬ 
formance  analysis  of  trackers  also  apply  here. 

I 
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2.6  USE  OF  THE  ANALYTIC  TECHNIQUES 

The  overall  approach  to  using  the  techniques  discussed  in  the  previous 
sections  is  summarized  in  Figure  2-3.  The  first  step  is  to  define  the  particu¬ 
lar  tracker  problem  with  the  result  being  a  set  of  system  models,  including  a 
dynamics  model  and  a  measurement  model.  Then  the  performance  bounds  and  approxi¬ 
mation  techniques  are  used  to  delimit  the  parameter  space.  Then  for  that  para¬ 
meter  space,  Kushner's  equation  can  be  solved  numerically  to  obtain  an  exact 
solution  of  the  bounds  for  comparison,  adjusting  the  delimited  parameter  space 
as  necessary  to  achieve  the  desired  results.  Then  algorithms  can  be  defined 
that  approximate  the  solution  to  Kushner's  equation,  representing  perhaps  the 
actual  tracker  implementation.  The  algorithms  can  then  be  compared  to  the 
exact  Kushner  solution  to  indicate  how  close  to  optimal  they  are.  The  proce¬ 
dure  is  iterated  until  acceptable  results  are  found. 
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FIGURE  2-3.  USE  OF  THE  ANALYTIC  TECHNIQUES 


SECTION  3.  TARGET  DYNAMICS  MODEL 


In  this  section,  a  model  is  developed  to  describe  target  and  target 
dynamics  for  a  three-dimensional  target  moving  with  6  degrees  of  freedom 
(3  translational,  3  rotational).  This  model  is  specialized  to  the  primary 
case  of  interest  in  this  study:  A  tactical  vehicle  (tank)  where  motion  is 
constrained  by  local  terrain  and  maneuvers  are  constrained  by  available 
controls  through  interaction  with  the  local  terrain.  The  target  contributes 
to  the  state  vector  and  dynamics  stochastic  equation.  A  simplified  target 
model,  approximated  by  a  rectangular  solid,  is  defined  and  the  basic  features 
of  the  resulting  two-dimensional  image  projected  into  the  image  plane  are 
discussed  in  Appendix  B.  This  model,  both  target  and  dynamics,  when  combined 
with  the  scene  model  developed  in  section  4,  provides  the  stochastic  differential 
equation  model  for  the  tracker  estimation  problem.  Also  presented  in  Appendix  D 
are  typical  numerical  values  for  the  parameters  that  determine  the  basic 
maneuvering  capability  of  tanks. 

3.1  THREE-DIMENSIONAL  TARGET  MODEL 

The  geometry  for  the  three-dimensional  target  model  as  related  to  the 
sensor  image  plane  is  shown  in  figure  3-1  for  a  rectangular  solid  target. 

The  concept  involved  is,  however,  applicable  to  arbitrary  three-dimensional 
target  shapes.  The  basic  approach  is  to  define  the  orientation  of  the  target 
in  terms  of  the  transformation  from  target  body  fixed  coordinates  to  sensor 
coordinates.  The  three  angles  necessary  to  define  this  transformation,  along 
with  the  target  position  relative  to  the  sensor,  provide  a  mapping  of  the 
three-dimensional  target  shape  to  a  two-dimensional  apparent  shape  via  sensor 
projection  model  developed  in  section  4. 


3-1 


The  target  body  coordinate  axes  x,  ,  y  ,  z,  are  defined  so  that  x  lies 

ODD  b 

along  the  direction  of  the  principal  force  used  for  speed  control,  z  along 

b 

the  principal  axis  for  turning  or  direction  changes,  and  yfa  forming  a  right- 
handed  coordinate  system.  The  two-dimensional  shape  is  determined  by  the 
aspect  angle  e  (yaw)  and  the  pitch  angle  P,  while  the  roll  angle  $  determines 
the  orientation  of  the  shape  in  the  image  plane.  The  transformation  from  body 
coordinates  to  sensor  local  vertical  coordinates  (x,y,z)  is  determined  by 
the  following  three  rotation  transformations. 

1).  Target  Yaw 

sin  6  O' 
cos  9  0 

0  1J 


H3- 


2) 


cos  9 
| -sin  0 
0 

Target  Pitch 
1 


Mi  - 


3). 


Image  Roll 

"cos  $ 


H 


0 

sin  4 


0 

cos  P 
-sin  P 

0 

1 

0 


0 

sin  P 
cos  P. 

-sin  $ 

0 

cos  $  J 


The  complete  transformation  is  then  given  by 

T-  K2  Si  Ms 


The  remaining  three  parameters  for  the  target  model  are  the  x^,,  yT>  z posi¬ 
tion  of  the  target  relative  to  the  sensor. 
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3.2  DYNAMICS  MODEL  FOR  A  TANK 

The  basic  model  for  a  maneuvering  tank  is  depicted  in  Figure  3-2.  The 
salient  features  of  this  model  are  quite  simple.  First,  the  tank  is  confined 
to  follow  the  local  terrain,  excluding  momentary  bump  effects.  Second,  the 
only  possible  maneuvers  are  linear  acceleration  along  the  hull  major  axis 
and  a  turn  around  what  is  basically  the  turret  axis,  and  these  are  only 
possible  through  interaction  with  the  local  terrain.  Thus,  the  possible 
direction  of  a  maneuver  depends  on  the  tank  orientation  and  the  local  terrain 
slope.  In  addition,  tank  type  restricts  the  possible  maneuvers,  e.g.,  maximum 
turn  rate,  speed,  acceleration.  Knowledge  of  maneuver  capability  is  essential 
for  obtaining  a  good  dynamic  model  for  filtering  and  prediction.  Similar 
techniques  are  used  for  other  vehicle  types. 

The  basic  equations  of  motion  for  this  tank  can  be  written  in  the  form 


where 

F  is  motor  or  braking  force 

Kj  and  K2  are  damping  constants 

t  is  the  torque  force  causing  vehicle  turn 

’  is  the  turning  rate 
and 

m  and  I  are  the  tank  mass  and  moment  of  inertia  about  the  turning 
axis. 
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FIGURE  3-2.  SAMPLE  TARGET  MODEL 


From  the  definitions  of  the  transformation  from  target  to  sensor  coor¬ 


dinates  above,  the  components  of  the  vector  j  can  be  defined  as 

Cos  <j>  Cos  o  -  Sin<j>Sin  o  Sin  P 
J  -  Sin  e  Cos  <p 

Cos  8  Sin  p  +  Cos  <p  Sin  0  Sin  P 

when  expressed  in  sensor  local  vertical  coordinates. 


There  are  two  types  of  disturbances  associated  with  this  dynamic  model. 

The  first  type  termed  controlled  disturbances  are  the  result  of  the  maneuvering 
by  altering  the  force  F  and  torque.  The  second  type  is  the  uncontrolled 
effects  induced  by  terrain  variations,  which  are  primarily  small  disturbances 
on  the  pitch  and  roll  of  the  target.  The  model  for  the  torque  disturbance 
is  shown  in  Figure  3-3,  and  consists  of  two  types:  Generalized  Poisson 
which  models  random  turning  of  the  tank  under  operator  control  and  white 
Gaussian  representing  terrain  induced  effects. 

The  target  portion  of  the  dynamics  state  equations  then  can  be  written 
as  follows: 


j-6 


X 


*o  =  y 


POSITION 


x3  =  z 


x4  =  X 


Xc  =  y 


VELOCITY 


X6  =  Z 


x7  =  F/m 


DRIVING  FORCE 


xg  =  P 


X10  =  9 


X11  =  9 


X12  T 


ANGLES 


ANGLE  RATE 


Torque 


Using  these  definition  along  with  the  equations  of  motion  and  the 
disturbances  gives  the  following  set  of  state  equations 
dx3  =  x^dt 
dx^  =  x^dt 
dx3  =  x6dt 

dx4  =  ^(cos  Xg  cos  Xj0  -  Sin  Xg  Sin  X^Q  Sin  Xg)  -  ^1_  X^j  dt 

dx5  =  [X7(-Sin  X10  cos  Xg)  -  ^1  Xg  j  dt 

m  J 

dx6  =  CX7^cos  X1Q  Sin  X8  +  C0S  X8  Sin  X10  Sin  X9^  '  —  Xgl  dt 


d x7  -  dF^  +  “  F-j  Wiener,  Fg  Poisson  processes 

dXg  =  d$j  -  Wiener  driving  process  on  roll 

dXg  *  dP  -  Wiener  driving  force  on  pitch 


This  general  form  for  the  target  portion  of  the  dynamic  equation  is 
extended  with  the  case  of  the  clutter  and  measurement  models  in  the  next 
section  to  provide  the  complete  stochastic  differential  equation  formulation 
of  the  target  tracking  problem. 


SECTION  4.  SCENE  MODEL 


In  order  to  get  useful  bounds  on  tracker  performance,  a  scene  model 
which  includes  the  effects  of  projecting  a  three-dimensional  scene  onto  a 
two-dimensional  plane  is  essential.  The  model  must  also  include  the  effects 
of  target/clutter  occlusion,  blurring,  spatial  sampling,  and  noise  degrada¬ 
tions.  Any  scene  model  which  does  not  include  these  effects  is  unlikely  to 
give  useful  results. 


4.1  PROJECTIONS 

Projection  effects  are  most  easily  handled  by  using  the  method  of 
perspective  transformations  to  transform  a  given  point  on  the  object  of 
interest  to  the  image  plane.  The  transformation  equations  are  contained  in 
any  book  on  scene  analysis  or  computer  graphics  (for  example,  ref.  [39]). 
Application  of  these  methods  to  our  problem  follows: 

Let  r  =  [x,y,z]T  be  the  vector  coordinates  of  a  point  in  space  relative 
to  a  sensor  set  of  fixed  coordinates  (see  figure  4.1).  The  sensor  line  of 
sight  axis  is  assumed  to  lie  along  the  z  axis  (third  coordinate).  Then  if 
the  sensor  has  focal  length  f,  and  f  '<  z: 


(4-1) 


where  u,v  define  (in  focal  length  units)  the  projection  of  the  point  on  the 
image  plane.  Now,  if  this  point  is  represented  as  F'  ®  Qx'.y'.z'J  in  some 
other  fixed  coordinate  frame  rotated  and  translated  from  the  sensor's 


(4-2) 


where  rQ  represents  the  offset  between  the  two  coordinate  systems  and  T 
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FIGURE  4-1.  PROJECTION  GEOMETRY 


the  rotation  matrix  is  a  function  of  0y*0p»6r»  the  yaw,  pitch  and  roll  angles 
between  the  two  coordinate  systems.  By  substituting  (4-2)  into  (4-1),  the 
projected  image  coordinates  of  a  given  point  in  space  relative  to  an  arbitrary 
fixed  coordinate  system  can  be  computed.  If  an  object  in  physical  space  has 
its  surface  defined  by  the  equation 

q'(r’)  =  0  (4-3) 

and  on  this  surface  (assumed  to  be  perfrectly  Lambertian)  it  has  an  emitted 
I.R.  intensity  profile  I'(r'),  then  to  predict  the  sensed  image  we  can  again 
use  (4-1)  and  (4-2)  to  produce: 

q(r)  =  q'(T-1(r-ro))  (4-4) 

I  (r)  =  I'(T-1(r-r  ))  (4-5) 

^  O 

where  q(r)  and  I(r)  are  the  transformed  surface  and  intensity  functions, 
respectively.  Now,,  representing  q(r)  as  q(x,y,z),  we  can  produce  the  function 
Qz(u,v)  defined  by 

!min  z  j  q(uz ,vz , z)  =  0|  if  3z  q(uz,vz,z)  =  0 

(4-6) 

-1  if~3z  q(uz,vz,z)  =  0 

The  interpretation  of  Qz(u,v)  is  that  it  represents  the  z  axis  coordinate  of 
the  "closest"  part  of  the  object  at  any  given  set  of  projection  coordinates 
subject  to  the  approximation  of  (4-1).  It  is  negative  for  a  (u,v)  position 
that  does  not  image  the  target  (it  is  assumed  that  q(x,y,z)  *»  0  -*■  z  >  0) . 

Using  (4-6)  and  (4-1),  we  can  solve  for  the  idealized  blur  and  noise-free 
image  of  the  object,  J(u,v)  in  terms  of  its  transformed  intensity  profile 
I(x,y,z) . 


*' 

f 

.4 

* 
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J(u,v)  =  S(Q)  •  I(Q-u,Q-v,Q) 


(4-7) 


fo,  X<0 

where 

S(x)  = 

1  1,  XiO 

and 

Q  =  Qz(u,v)  . 

(4-8) 


In  practice,  J(u,v)  is  a  difficult  function  to  compute  analytically. 

Appendix  B  treats  a  related  set  of  transformations  for  a  rectangular  solid. 

Another  example  could  be  an  ellipsoid  as  defined  in  Table  4-1.  An  ellipse 

is  a  good  general  shape  and  as  shown  in  the  table,  can  be  represented  in 

body  coordinates  by  elliptical  equations  defined  by  the  three  principal  axis 

radii,  r  ,  r  ,  r  .  If  the  temperature  variation  over  the  ellipsoid  is 
x  y  z 

negligible,  the  intensity  over  the  object  is  constant  and  the  transformation 
of  the  object  to  line-of-sight  coordinates  can  be  accomplished  as  indicated 
in  Table  4-2,  where  rQ  is  the  distance  from  the  object  to  the  LOS  coordinates 
located  at  the  sensor  and  T  is  the  yaw-pitch-roll  transformation  defining  the 
orientation  of  the  object.  Since  the  object  has  uniform  intensity,  the  3- 
dimensional  "intensity"  in  LOS  coordinates  is  just  the  transformation  of  each 
point  in  the  object.  The  sensed  2-dimensional  image  can  be  found  by  use  of  a 
"nearest"  function  as  also  shown  in  Table  4-2.  This  function  provides  a  means 
of  formalizing  the  concept  that  the  front  of  the  object  blocks  out  all  points 
in  the  object  behind  it.  The  intensity  function  J  given  in  the  table  then  is 
the  observed  object  since  it  is  1  for  front  points  and  zero  for  points  off 
the  object.  This  concept  is  extended  to  include  blocking  by  clutter  in  the 
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next  section. 


TABLE  4-1.  EXAMPLE  OF  AN  OBJECT  MODEL 


WHERE  r'  =  [x'. 


4.2  TARGET-CLUTTER  OCCLUSIONS 

The  effects  of  target  and  clutter  occlusions  can  now  be  included.  Let 
Qzt  be  the  Q  function  (reference  eqn.  (4-6))  for  the  target,  and  Qzc(u,v) 
represent  the  equivalent  Q  function  for  the  object  composed  of  the  aggregation 
of  all  potentially  obscuring  clutter  objects.  Similarly,  let  JT(u,v)  and 
Jc(u,v)  represent  the  unobscured  target  and  clutter  intensity  projections, 
respectively,  as  in  (4—7) . 

Define  the  function  D(u,v)  such  that 

II  if  QZ(J(u,v)  <  0 

0  if  QZT(u,v)  <  0  (4-9) 

S(Qzc(u,v)  -  QZT(u,v))  otherwise 

with  S(x)  as  defined  in  (4-8).  Loosely  speaking,  D(u,v)  *  1  when  the  target 
obscures  the  clutter  at  image  position  (u,v)  and  0  otherwise. 

The  observed  idealized  (blur  and  noise-free)  intensity  of  a  target 
moving  in  and  out  of  clutter  is  then  given  by: 

JTC^U,V0  ■  D(u,v)  •  JT(u,v)  +  ( 1-D(u ,v) )  Jc(u,v)  (4-10) 

as  can  be  seen  in  Figure  4-2. 
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A. 3  BLUR  AND  SAMPLING 


First  order  blurring  effects  are  included  by  assuming  a  uniform  2-D  blur 
function.  In  general 


ho(u’V>  =  ~2  f  f  JTC(a>6) 

£  —  ao  _ ao 


b^(u-a,v-B)  dadB 


(4-11) 


where  the  blur  function  b^x.y)  is  given  by  the  system  PSF.  In  this  report 
we  will  commonly  use: 

(1  if  |x|<A/2, |yjsf 

bA(x,y)  =  rect(x/A  ,y/A)  =  <  (A-12) 

(0  otherwise 

hQ(u,v)  can  then  be  sampled  at  points  (u^v^)  to  produce  the  vector  observation 
function  h 


h=  hQ(u2,v2),  ho(u3,v3)...]T  (4-13) 

Sampling  intervals  can,  of  course,  vary.  In  many  examples  in  this  report  we 
typically  use  one  sample  per  blur  element  in  a  uniform  rectangular  grid 


u±  =  m(i) 
v,  *  n(i) 


ro(*)»n(")  integer  functions 


( 4— 14 ) 


A. A  MEASUREMENT  MODEL 


Once  we  have  selected  the  h  function,  which  has  as  its  implicit  variables 
the  state  vector  describing  the  target  and  clutter  objects  being  imaged,  noise 
is  added  to  the  model  to  represent  sensor  noise  in  the  observation.  We  will 
model  the  noise  to  be  Gaussian,  and  temporally  and  spatially  white.  This  is 
represented  by 

dy  =  h(x)dt  +  dz  (A— 15) 


where  x  represents  the  augmented  state  vector  which  includes  vehicle  dynamics, 

clutter  states  and/or  obscuration  variables  (states  for  producing  D(u,v)  of 

eqn.  (4-9)),  and  y  is  the  observation  output  for  this  model. 
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The  covariance  of  the  noise  is  related  to  the  noise  effective  temperature 


of  the  sensor  and  frame  rate  by: 

E[dz  d2T]  -  I  it  ■  (4-16) 

where  E  Clis  expectation, 

NEAT  is  the  sensor  noise  effective  temperature, 

B  is  the  frame  rate 
I_  is  the  identity  matrix. 

Note  that  the  output  of  an  actual  sensor  would  be 

t 

I(t)  =  f  dy  =  y(t)  -  y(t-A)  (4-17) 

t-l/B 

and  is  not  to  be  confused  with  I(r)  of  equation  (4-5).  While  the  equa¬ 
tions  given  by  (4-17)  are  the  equations  actually  satisfied  by  any  real  system, 
the  equations  of  (4-15)  and  (4-16)  are  more  convenient  for  mathematical  analysis 
and  contain  exactly  the  same  statistical  information*.  Hence,  any  results  ob¬ 
tained  when  the  measurements  are  specified  by  (4-15)-(4-16)  apply  when  they 
are  specified  by  (4-17). 

Figure  4.3  illustrates  the  scene  modeling  process  described  in  this 
section. 


*In  the  limit  as  &+<»  (fast  frame  rate). 
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FIGURE  4-3.  SCENE  MODELING 


5.0  ANALYTIC  RESULTS 


In  this  section,  we  introduce  two  analytical  methods  for  evaluating 
electro-optical  image  trackers,  namely,  linearization,  and  Bobrovsky-Zakai . 

Both  methods  consist  of  replacing  the  non-linear  tracker  model  with  an  equivalent 
linear  model,  from  which  a  matrix  Riccati  equation  is  derived.  A  bound  or  an 
approximation  to  the  mean  square  tracking  error  is  then  obtained  by  solving 
the  matrix  Riccati  equation. 

In  section  5.1,  the  Bobrovsky-Zakai  method  is  used  to  derive  a  lower 
bound  to  the  RMS  tracking  error  for  a  simple  rectangular  target  constrained 
to  move  horizontally  across  the  FOV.  At  any  given  time,  the  target  could  be 
occluded  by  a  stationary  clutter  object.  The  dynamical  model  used  is  simple. 

In  section  5.2  and  5.3,  the  linearization  method  is  used  to  approximate 
the  RMS  tracking  error  for  more  realistic  dynamics  models. 


5.1  BOBROVSKY-ZAKAI  BOUND 

In  principle,  the  Bobrovsky-Zakai  method  could  be  applied  to  the  general 
dynamics-scene  model  introduced  in  sections  3  and  4.  In  section  3,  the  general 
dynamics  model  was  given  as 

dtf  =  f(x,t)  dt  +  C  dw,  (5-1) 

where  3?  is  a  11-dimensional  state  vector  specifing  the  target's  location, 

velocity,  yaw,  pitch,  roll,  and  angle  rates. 

The  general  scene  model  given  in  section  4  conveniently  specifies  the 
measurement  equation  as 

dy  -  (xA)  dt  +  D  dz  (5-2) 

where  x 

A  represents  the  augmented  state  vector  which  includes  vehicle 
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dynamics  (the  t  vector  in  equation  5.1)  ,  and  clutter  states  and/or  obscuration 
variables,  y  is  the  observation  vector.  The  notation  used  here  is  slightly 
different  than  that  used  in  section  3  and  4.  The  notation  change  is  only  to 
facilitate  the  application  of  the  Bobrovsky-Zakai  method.  Matrices  C  and  D 
are  related  to  the  forcing  function  and  measurement  noise  of  sections  3  and  4 
as  follows: 

CCTdt  =  E  (dw  dST}  (5-3) 

where  w  is  the  random  forcing  function  in  section  3,  and 

DDTdt  =  E  {dz  dzT}  (5-4) 


where  t  is  the  measurement  noise  in  section  (4). 

For  the  general  tracking  problem  specified  by  equation  (5-1)  and  (5-2), 
the  Bobrovsky-Zakai  bound  can  be  obtained  by: 


(1)  Set 


(5-5) 


where 


3f 

3X 


is  the  vector  gradient  of  f 


I 


Solve 

for 

B  in 

bW)'1 

B  =  e{( 

(4-*)T«*Tri  • 

(4- a)! 

(5-6) 

3x 

3x 

+  E  ' 

3*T 

(ddt)‘ 

-1  3h  ] 

’ 

1 

1  3X 

l 

3x 

, 

3h 

where  gi  ’s  vector  gradient  of  h 

(3)  Replace  the  nonlinear  problem  with 

dGT  =  A  "u  dt  +  C  dw  (5-7) 

did  =  B  u  dt  +  D  d?  (5-8) 
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(4)  Solve  the  matrix  Riccati  equation  for  the  linear  system  specified 
by  equation  (5-7)  and  (C-8),  which  is  given  by 

dV  =  AV  +  VAT  +  CCT  -  VBTBV  (5-9) 

dt 

The  mean  square  error  V  given  by  the  solution  to  the  Riccati  equation  will 
then  be  a  lower  bound  to  the  mean  square  error  of  the  nonlinear  system  (5-1),  (5-2). 

In  general,  equations  (5-6)  and  (5-7)  cannot  be  evaluated  analytically  for 
a  realistic  dynamics-scene  model;  however,  both  can  be  easily  evaluated  by  Monte- 
Carlo  methods.  The  computational  burden  involved  in  simulating  equations  (5-6)  and 
(5-7)  will  in  general  be  less  than  that  required  to  simulate  the  optimal  tracker 
via  Kushner's  equations  or  the  integral  method  described  in  section  6.1. 

It  has  been  pointed  out  that  for  any  realistic  problem,  an  easily 
implemented  Monte  Carlo  simulation  will,  in  general,  be  required  to  evaluate 
certain  averages  which  are  used  in  the  Bobrovsky-Zakai  bound.  However, 
frequently  a  qualitative  understanding  of  tracker  performance  versus  system 
parameters  can  be  obtained  by  analyzing  simple  scenarios  for  which  analytic 
bounds  can  be  obtained. 


Hence,  the  general  dynamics-scene  model  will  be  restricted  to  the  simple 
scenario  where  a  rectangular  target  with  fixed  size  and  aspect  ratio  moves 
horizontally  across  the  FOV  with  a  Gaussian  white-noise  velocity.  It  is  further 
assumed  that  there  is  a  stationary  clutter  object  in  the  FOV  which  at  any 
time  might  partly  or  completely  occlude  the  target.  (See  figure  5.1) 

For  this  simple  scenario,  equation  (5-1)  reduces  to 


d  Xy  =/q  dw 

(target  motion) 

(5-10) 

“"c-o 

(clutter  location 

is  fixed) 
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FIGURE  5.1.  SCENE  MODEL 


where  q  is  the  target's  velocity  spectrum  amplitude.  The  measurement  equation 
(5  -2)  reduces  to 

dy  =  (x j ,  xc)  dt  +  D  dz  (5-11) 

where  for  pixel  i 

hi(xT,xc)  =  Tj  hT  (xT)  (1-  hc  (xc))  +  Tc  hc  (xc) 


=  target  intensity  offset 
Tc  =  clutter  object  intensity  offset 

hT  (xT)  =  target  indicator  function 

h  (x  )  =  clutter  indicator  function 
C  •  c 


The  app’ication  of  the  Bobrovsky-Zakai  method  results  in  a  RMS  error  bound 
given  by 


RMS  lower  bound  = 


rq/aT ( 1 -  ex  p ( - /aTq/r' t)) 
(1  +  exp(-/aj(l/rl 


H 


(5-12) 


where 


qco  =  density  of  clutter  pixels  =  #  pixels  on  clutter  object/ 
#  pixels  in  FOV 


qce  =  density  of  vertical  clutter  edge  pixels 
clutter  edge/#  pixels  in  FOV 


(NEAT)* 


B  =  frame  rate 


=  #  pixels  on  vertical 
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Ty=  target  temperature  offset, 
q  =  specturm  amplitude  of  velocity  process 
A  =  blur  size 


Np  =  number  of  blur  elements  across  target  in  the  vertical  direction. 
A  detailed  derivation  of  this  result  is  given  in  appendix  E. 

Figure  5.2  demonstrates  the  variation  of  this  bound  as  a  function  of 
fraction  of  effective  FOV  covered  by  the  clutter  object,  where  effective 


FOV  is  defined  as  that  part  of  the  FOV  for  which  it  is  possible  for  the 


target  to  lie. 
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ERROR  (PIXELS) 


EFFECTS  OF  CLUTTER  ON  IMAGE  TRACKING 
FIGURE  5.2 


0  SNR  =  AT/NEAT  =  1 
o  B  *  30HZ 

oq«3  (PIXELS/SEC. )2/HZ 
o  TARGET  SIZE  =  5  x  5  PIXELS 

o  TARGET  AND  CLUTTER  POSITION  ASSUMED  KNOWN  AT  TIME  0. 


'•2.  r 


I 

■6 

oo 

S  >2 


75%  CLUTTER 


CLUTTER 
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- - ,  .  ■  f - !  ■  r  . 7  "  "t  r  i  t 

*2  *5  <4  .S  .7  >8  /. 
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5.2  APPLICATION  OF  LINEARIZATION  TO  I-D  FORCE  MODEL 


A  second  approach  to  performance  analysis  of  imaging  trackers  is  to 
find  approximate  results  by  linearization.  Starting  with  the  simplest 
dynamics  model  these  linear  equations  are  solved  to  provide  analytic 
relationships  between  the  RMS  tracker  error  and  the  appropriate  descriptive 
parameters.  The  results  are  extended  to  progressively  more  general  dynamics 
models  to  find  generalizations  which  relate  desired  figures  of  merit 
(e.g.  RMS  error)  to  system  parameters. 

Use  of  linearization  results  in  an  extended  Kalman  filter,  which  can 
then  be  solved  to  get  an  approximate  solution.  This  section  shows  how  that 
approach  can  lead  to  useful  functional  relationships  and  point  the  way  to 
generalizations  of  particular  interest  in  the  error  covariance  matrix  V, 
which,  after  linearization,  satisfies  a  Riccati  equation  of  the  form: 

^  «  AV  +  V  AT  +  M  -VNV  (5-13) 

at 

where  A  is  the  dynamics  matrix,  M  is  a  stochastic  disturbance  matrix  and  N  is 
a  measurement  matrix. 

Three  cases  of  the  equation  for  a  linearized  Kalman  filter  with  one¬ 
dimensional  motion  have  been  analyzed,  and  closed-form  solutions  have  been 
obtained  for  each  case.  In  each  of  these  cases,  the  dynamics  equation  becomes 

linear,  so  that  the  linearization  process  need  be  applied  only  to  the  measurement 
equation. 

The  one-by-one  case  is  the  Brownian  motion  and  the  target  position  is 

regarded  as  a  one-dimensional  process. 

The  second  case  is  a  two-state  problem  (two-by-two  matrices)  where  the 
acceleration  is  white  noise,  so  that  the  motion  is  affected  by  an  integrated 
process. 

The  third  case  is  a  three-state  problem  (three-by-three  matrices)  where 
the  jerk  (i.e.,  the  time  derivative  of  acceleration)  is  a  white  noise  process. 
Thus,  the  motion  is  affected  by  a  doubly  integrated  process. 

The  velocity,  acceleration,  and  jerk  spectral  densities  for  the  three 
cases,  respectively,  Q^,  Q^.  and  Q^,  are  defined  in  Table  3-2,  and  define  the 
M  matrix  for  each  case.  Note  how  the  time  dependence  changes  for  each  case, 
and  that  the  reciprocal  of  the  observation  jitter  spectral  density,  q  , 
appears  as  the  upper  left  element  in  each  of  the  N  matrices. 
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The  jitter  spectral  density,  q,  may  be  shown  to  be,  in  effect,  the  same 

quantity  developed  in  the  ATAC  Autocuer  Modeling  Study  report  [30]  using  the 

Cramer-Rao  bound  to  estimate  the  uncertainty  in  estimating  the  position  of  a 

target.  This  can  be  seen  as  follows:  The  quantity  q  is  a  product  of  three 

2 

factors,  of  which  the  first,  (6/SNR)  ,  is  the  variance  of  a  single  edge  pixel. 

The  second  factor ,  6/W,  the  reciprocal  of  the  number  of  pixels  on  a  target 

edge,  gives  the  spatial  integration  effect  in  reducing  the  edge  variance. 

The  third  factor,  1/B,  relates  the  jitter  amplitude  to  a  jitter  spectral 
density. 

The  solutions  to  the  equation  for  the  three  cases  are  given  in  Table  5-3. 
Note  that  for  the  steady-state  case,  the  observation  spectral  density  q  appears 
as  a  linear  factor  in  each  case  multiplied  by  the  reciprocal  of  a  certain 
time  constant.  In  fact,  the  steady  state  for  all  three  cases  is  given  by 
the  formula 

V, .  (“)  *  q  *  T  1  (5-14) 

ll  n 

where 

Tn  =  2_(n~1)'/2  (q/Qn)1/(2n) 
for  n  *  1,2,3,  respectively. 
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TABLE  5-1.  RICCATI  EQUATION  APPLIED  TO  VARIOUS  MODELS 


TABLE  5-2.  DEFINITIONS  OF  MATRICES  APPEARING  IN 
RICCATI  EQUATION  FOR  LINEARIZATION  EXAMPLES 


TABLE  5-3.  RESULTS  OF  LINEARIZATION  EXAMPLE 


Brownian  Motion  Variance: 
Vn(co)  =  q/T^1 


White  Noise  Force  Variance: 


White  Noise  Jerk  Variance: 


Details  of  the  development  of  the  solution  of  second  and  third  cases  are  given 
in  Appendices  C  and  D.  Here,  the  procedure  is  applied  to  the  first  case  to 
illustrate  this  method.  In  this  case,  the  equations  reduce  to  two  scalar 
equations : 


dx  =  dw 

dy  =  h  (x)  dt  +  dz 

where  dw  and  dz  are  Wiener  Processes. 


(5-15) 

(5-16) 


The  linearization  applied  to  the  measurement  is  the  scalar  degenerate 
version  of  the  linearization  given  in  Appendix  C  arJ  the  variance  V  obeys 
the  equation. 

~  =  M-VNV  =  Q]  -  V2  (5-17) 

Where  Qj  and  q  are  as  defined  in  Table  5-2.  The  solution  of  this 
scalar  differential  equation  from  an  initial  condition  V  =  0  is 

V  =  a  tanh  bt  rs-isi 


Direct  substitution  of  this  solution  into  the  differential  equation 
provides  a  solution  for  the  constants  a  and  b  as: 

a  =  v/  Q-j  q 

b  =  v/Q-j/q 

so  that 

V  =  /Ol?  tanh  /(^/q't  (5-21) 


(5-19) 

(5-20) 
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5.3  APPLICATION  OF  LINEARIZATION  TO  TWO-DIMENSIONAL  FORCE  MODEL 

This  model  is  like  the  one-dimensional  model,  except  that  the  target  is 
permitted  to  move  in  two  perpendicular  directions  with  random  white- 
noise  force  components  in  both  directions.  This  is  therefore  a  4-state 
problem  in  which  the  state  variables  are  defined  in  terms  of  the  x  and  y 
components  and  velocities  as  follows: 

x^  =  x 

X2  =  X 

x3  =  y  (5-22) 

\  =  y 

Thus  the  covariance  matrix  has  sixteen  (equals  four  times  four)  components,  of 
which  six  are  the  same.  Thus  one  is  led  to  a  Riccati  equation  equivalent  to 
ten  coupled  non-linear  differential  equations. 

When  one  imposes  the  requirement  that  the  covariance  matrix  elements  are  all 
to  be  zero  initially  at  time  zero,  one  finds  upon  examining  the  equations 
that  the  matrix  elements  correlating  x  and  y  positions  and  velocities,  namely 
V^3  V23  and  will  remain  zero  for  all  time,  since  they  have  no 

non-zero  forcing  terms.  Thus  the  x  and  y  equations  decouple  into  two  sets  of 
three  equations  each.  Each  set  is  equivalent  to  the  set  of  coupled  non-linear 
differential  equations  solved  in  the  preceding  section.  Thus  the  solutions 
are  essentially  the  same,  except  for  x  and  y  parameters  which  may  be  different 
for  the  two  directions. 

For  example,  the  steady  state  tracking  error  for  the  x  and  y  directions 


is  given  by 


(5-23) 


°x  '  <4  «x'%3>'/8 
",  - 14 


(5-24) 


where  q  «  ^  (SNR)2  B  (5-25) 

X  6J 

as  before  with  W  the  target  width,  but 

q  «  ^  (SNR)2  B  (5-26) 

y  r 

with  L  the  target  length.  Even  if  Q  and  Q  are  equal,  the  tracking 

x  y 

error  for  the  two  directions  will  be  different  if  the  target's  length  and  width 
are  unequal. 
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SECTION  6.  SIMULATION 


6.1  Introduction 

The  purpose  of  the  simulation  is  two-fold.  First,  it  provides  direct 
verification  of  the  validity  of  some  of  the  calculations  performed  in  the 
preceding  section.  Secondly,  it  provides  a  means  for  directly  extracting 
characteristics  of  the  probability  density  of  more  complicated  systems, 
whose  performance  may  be  more  difficult  to  estimate  in  a  tractable  fashion. 


6 . 2  Approach 


The  approach  that  was  employed  in  the  computer  simulation  is  to  first 
approximate  the  continuous  time  differential  models  with  appropriate  discrete 
time  difference  models,  i.e., 


x(K)  -  x(K-l)  =  ?(?(K-1))  +^(<-1) 

yOO  -7(k-D  »  7oT(km))  +'w2(k-i) 


or  equivalently  letting 

ZOO  =  y(K)  -  y(K-l) 
7(x(k-D)  =  "s(x(k-i))  +  -k(k-D 

we  have  from  (6-1) 

TOO  *  f(T(K-l))  +T1(k-i) 
Too  =  “h(x(K-l ))  +"w2(K-1) 
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E 


where  x(k) ,z(k)  are  intended  to  represent  the  state  and  differenced  observation 
vectors  respectively  at  the  kth  time  epoch,  i.e.  t  =  kAt.  Wj(k)  and  w^Ck)  are 
multivariate  independent  (each  dimension  and  for  each  k)  Gaussian  random 
variables. 

Under  this  set  of  assumptions  J^(u),  the  probability  density  of  the 
th  ^ 

state  vector  at  the  k  interval,  given  observations  of  ^(i)  for  i  <_  k  is  given 
by  (ref  4): 


P-x(K)|T(“lT(i,>  1SK>  5  JK<"> 


JK(u)  = 


Uu) 


(6-3) 


/  lk(“) 


where  Lk(u)  (an  unnormalized  probability  density)  is  given  by 


LK(U)  =  N(7(K)  -  h(u).  R)  X 


where 


f N((u-?(v)),Q)  JK.-,(v)  dv 


R  =  E[W2(K)  W2(K)] 
Q  =  E[W1 (k)  wJ(K)] 


(6-4) 


(6-5) 


are  the  noise  covariance  matrices  and  N(A,B)  is  proportional  to  the  Gaussian 


multivariate  density 


N{I,B)  =  exp  1  (AT  B'1  A)] 


(6-7) 


Figure  6.1  illustrates  the  operation  of  the  simulation.  In  the  first 
stage,  prior  to  actual  execution  of  the  simulation,  subroutines  representing 
the  "f  and  h  functions  are  linked  in  the  remainder  of  the  simulation  to  produce 


t 


an  executable  object  file.  Following  this  stage,  execution  begins  with  the 
routine  parameters  being  entered  from  the  keyboard.  These  include  output 
modes,  number  of  iterations  (N) ,  dimensionality  of  state  and  observation 
vectors,  Q  *  diagonal  elements  (since  these  are  assumed  diagonal  matrices), 
number  of  samples  of  "u  per  dimension,  starting  state  for  the  system,  and  initial 
state  probability  density  parameters. 

After  this  Jo(u),  which  is  represented  by  values  on  a  uniform  sampling 
grid,  is  initialized  as  is  the  initial  system  state. 

Then  for  each  iteration  of  the  simulation,  the  system  producing  observa- 
tions  is  simulated  via  the  f  and  h  subroutines  and  a  random  number  generator 
to  produce  simulated  noise.  Following  this,  the  next  L^Cu)  function  is  again 
computed,  point  by  point,  on  the  same  sampling  grid  that  Jkl(u)  was  defined 
on,  using  numerical  integration  techniques.  It  is  in  this  stage  that  most  of 
the  simulation's  throughput  is  used.  After  L^(u)  is  computed,  it  is  normalized 
to  produce  the  new  Jk(u)  and  the  results  are  printed  out  on  the  terminal  and 
a  file.  The  sampling  grid  is  then  redefined  to  accommodate  the  new  probability 
density  (allowing  for  spreading  and  shifts)  and  the  next  iteration  begins. 

Following  the  execution  stage,  a  graphing  program  can  read  the  files  con¬ 
taining  the  results  and  produce  graphs  of  time  tracks  of  estimated  means, 
standard  deviations,  or  probability  densities  to  allow  more  readily  interpret¬ 
able  results. 


6-4 


6.3  Simulation  Comparison 


Due  to  the  stochastic  nature  of  the  simulation,  any  quantity  observed 

during  a  single  run  of  the  simulation  should  be  considered  as  the  outcome  of 

a  random  experiment  and  combined  with  the  results  of  other  simulation  runs 

to  obtain  an  estimate  of  the  quantity  of  interest.  Here  the  quantity  of 

interest  is  RMS  tracking  error,  and  it  is  estimated  from  the  output  of  the 

simulation  as  follows: 

/n,  /  N 

RMS(t.)  -  52 

1  \k-l 

where  W  “  W  '  W’  =  N  ^  ek(ti)  ’ 

k=l 

N  is  the  number  of  simulation  runs,  the  target's  location  at  time 

t^  for  the  kth  simulation  run,  and  ^(t^  is  the  corresponding  optimal  esti¬ 
mate  of  X^t^). 

In  figure  6.2,  both  RMS(t^)  and  the  Bobrovsky-Zakai  bound  are  given  for 
2 

a  5  x  5  pixel  target  which  is  moving  horizontally  across  the  FOV  with  a 

Gaussian  white  noise  velocity.  Although  the  target  is  square  these  results 
hold  equally  well  for  rectangular  targets  since  the  motion  is  one  dimensional. 

The  height  of  the  target  is  just  used  to  build  up  SNR.  The  target's  velocity 
spectrum  amplitude  is  assumed  to  be  3  (pixel/sec)  /Hz  and  SNR  =  1.  Also,  the 
«T  of  equation  (5-12)  would  be  one,  since  there  are  no  clutter  objects  in  the 
scene . 

In  figure  6.3,  RMS(t^)  and  the  RMS  error  obtained  by  the  linear  approxi- 

2 

mation  method  are  given  for  a  2  x  2  pixel  target  which  moves  horizontally  across 
the  FOV,  satisfying  the  dynamics  equation 
dX  =  V  dt 

dV  “YT  dw 

where  X  is  the  target  location  and  V  is  its  velocity.  SNR  is  assumed  to  be  1 

2  2 

and  the  spectrum  amplitude  of  the  target's  acceleration  is  1.9  (pixels/sec  )  /Hz. 
For  this  simple  example,  the  results  from  the  linear  approximation  agree  well 

with  the  simulation. 
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SNR  =  AT /NEAT  =  1 
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FIGURE  6-2.  BOBROVSKY-ZAKAI  BOUND 


TARGET  SIZE 


2x2  PIXELS 


SNR  =  AT /NEAT  =  1 
TARGET  ACCELERATION  SPECTRUM 
AMPLITUDE  =  1.9  (PIXELS/SEC") "/HZ 

B  =  30  HZ 


.2.2, 


LINEARIZATION 


SIMULATION 


0  0.40  o'.  50  0.80  1.00 

TIME  (SEC) 

FIGURE  6-3.  LINEARIZATION  METHOD 


In  order  to  run  the  discrete  simulation,  the  following  inputs  are 
required: 

1)  .  R  =  NEAT2 

2) .  Q  =  q/^f 

where,  for  the  B-Z  example  q  is  the  velocity  spectrum  amplitude  and  Af  is  the 
frame  rate.  For  the  linearization  example  q  is  the  acceleration  spectral 
amplitude.  The  values  of  all  parameters  necessary  to  run  the  simulation  for 
the  B-Z  example  and  the  linearization  example  are  given  in  tables  6-1  and 
6-2,  respectively. 


TABLE  6-1.  INPUTS  FOR  B-Z  EXAMPLE 


#  OF  TRIALS 

30 

TARGET  SIZE 

2 

5x5  pixels 

FOV 

25  x  5  Pixels^ 

DELTA  T 

1. 

NUMBER  OF  STATES 

1 

INVERSE  Q 

( 30) ~ 1 

INVERSE  R 

1 

INITIAL  TARGET  STATE 

ASSUMED  KNOWN 

NUMBER  OF  TIME  SAMPLES 

30 

TABLE  6-2.  INPUTS  FOR  LINEARIZATION  EXAMPLE 


#  OF  TRIALS 
TARGET  SIZE 
FOV 

DELTA-T 

NUMBER  OF  STATES 
INVERSE  Q 
INVERSE  R 

INITIAL  TARGET  STATE 
NUMBER  OF  TIME  SAMPLES 


2x2  pixels 
10  x  2  pixels^ 
1 
2 

(1.9/ 30) — 1 
1 

ASSUMED  KNOWN 
30 


SECTION  7.  APPLICATION  TO  PRECISION  AIMPOINT  MAINTENANCE 
The  analytical  methods  and  results  descirbed  in  this  report  have 
application  in  the  design  and  analysis  of  electro-optical  systems  requiring 
tracking  functions  to  be  performed.  A  key  applications  area  includes 
precision  aimpoint  maintenance  systems. 

Precision  aimpoint  maintenance  is  a  required  function  for  systems  such  as 
directed  energy  weapons.  In  this  application,  low  tracker  jitter  and  RMS 
error  from  aimpoint  is  a  requirement.  Since  small  deviations  from  the  target 
position  are  relevant  in  this  application  the  linearization  approach  outline 
in  section  5.3  will  be  used  in  the  following  example. 


7.1  Target  Model 

The  target  to  be  tracked  in  this  case  will  be  modeled  to  be  imaged  as 
a  rectangle  of  length  L  (radians),  width  W  (radians),  and  intensity  offset 
AT  (°C),  blurred  by  a  characteristic  blurring  function  of  dimension  r  (radians). 
In  other  words: 


where 


and  ®  is  convolution 


AT-rect(£p  .  ^)®rect(£,  £ 


rect  (x,y)  = 


1  . 

0  , 


(7-1) 

M  <  1/2,  iy I  <  1/2 

I x |  >  1/2  or  |y|  >  1/2 


I,  ,(x,y)  represents  the  perceived  intensity  of  a  target  located  at 
image  position^,  v)  detected  at  sensor  azimuth  x  and  elevation  y. 

The  motion  of  the  target  (relative  to  the  sensor  line  of  sight)  will  be 
modeled  by  second  order  linear  stochastic  differential  equations  in  the 
azimuth  and  elevation  directions 

du  =  u  dt 
du  =  d  w 

.  U  (7-2) 

dv  =  v  dt 

dv  =  d  wv 

wjth  w  w  being  uncorrelated  Wiener  processes, 
u ,  v 

7.2  Sensor  Model 

The  sensor  model  will  treat  the  sensor  as  providing  continuous  samples 
in  time  of  a  discretely  sampled  (in  space)  intensity  function.  The  sampling 


*  Using  radians  is  equivalent  to  using  focal  length  units  of  section  4 
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interval  is  r  radians,  with  h  (odd)  samples  subtending  the  azimuth  and 
elevation  FOV.  Each  sample  will  be  corrupted  by  an  independent  white  noise 
source,  i.e. 

d  J  =  $  dt  =  dZ  (7-3) 

where  "3"  is  a  vector  of  length  n2  representing  the  observed  intensity 
samples 

~l  is  a  vector  Wiener  process 

E[dZ  c izT]  =  Tn  I  dt 

where  Tn  is  (NE A T)2/  B,  B  =  frame  rate 

4>  is  a  vector  representing  the  n2  samples  on  the  image  plane 

V  =  C*ir  *2,1 . *ml *  ^12’  ^22  •**  (7‘4) 

where 

=  f(i'T/2,  j '  •  r/2) 
i'  =  i  -  (m-1 )/2,  j’  =  j  -  (m-1 )/2 

and  f  (x,y)  is  the  focal  plane  intensity  (°C)  at  azimuth  x,  elevation  y. 

Sensor  motion  will  be  accounted  for  by  equation  (7-2)  since  the  angular  motion 
of  the  target  is  defined  relative  to  the  sensor  line  of  sight. 

7.3  Combined  Model 

Combining  the  sensor  and  target  models 

du  =  u  dt  dv  =  dt 

du  =  d  w„  dv  =  d  w 

u  v 

dJ  =  *(u,v)  dt  +  dZ 

V“-"  ■  V.  (  [< -  r«.  [  J  -  i^l]  rn) 

■u>v(x.y)  ■  rect  tSj  ©  rect  (x/r,  y/r) 
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(7-5) 


(7-5)  Is  an  example  of  a  linear  system  -  non  linear  observation 
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7.4  Performance  Analysis 


\ 

I 


As  previously  stated,  this  application  is  concerned  with  a  low  jitter 
and  RMS  error  from  aimpoint  requirement.  Thus,  the  results  of  the  linearization 
method  introduced  in  section  5.3  will  be  used. 


The  steady  state  tracking  error  (from  section  5.3)' 


is 

“el  ■  "VeIi 

for 

the 

elevation  direction 

(7-6) 

and 

°AZ  ■  <4  WaI11'8 

for 

the 

azimuth  direction 

(7-7) 

where 

QaZ’  ^EL  =  sPectrum  amplitude  of  the  angular  acceleration 
due  to  jitter 

qEL  =  3}  (SNR)2b 


qAZ  “  %  (SNR)2  B 

B  =  frame  rate 
6  =  blur  size. 


In  order  to  preceed  further,  the  following  sensor  parameters  are  assumed 
6  =  blur  size  =  180  yrad 
B  =  frame  rate  =  30  HZ 

QAZ*  QEL  =  18.7  (mrad/sec2)2/HZ 
The  target  size  is  assumed  to  be 

W  =  .66  mrad 
1=1.9  mrad 

The  elevation  and  azimuth  steady  state  RMS  tracking  error  for  these 
parameters  are  given  in  figure  7.1-  This  figure  could  then  be  used  to  set 
performance  requirements.  For  example.  If  an  aimpoint  accuracy  of  30p  rad 
Is  required,  then  It  Is  seen  that  SNR  »  AT/NE  AT  must  be  no  lower  than  2.5. 


I 
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SECTION  8.  CONCLUSIONS 


This  section  will  summarize  the  work  performed  under  this  contract  and 
draw  certain  conclusions  related  to  this  work. 

The  first  conclusion  one  can  draw  based  on  the  study  of  both  the  target 
and  the  scene  model,  as  well  as  the  definition  of  the  tracker,  is  that  the  tech¬ 
nique  of  applying  non-linear  stochastic  differential  equations  has  the  power 
to  completely  model  electro-optical  tracking  systems.  A  key  point  to  be 
noted  here  is  that  both  the  target  dynamics  as  well  as  sensor  noise  can  be 
described  using  the  same  types  of  Martingale  processes  and  can  be  treated 
in  a  mathematically  unified  manner.  The  model  that  has  been  developed  has  the 
ability  to  include  maneuvers,  clutter,  target  obscuration,  and  can  include 
sensor  blur,  noise,  and  sampling  effects  as  well.  This  ability  to  treat 
the  general  problem  in  a  unified  way,  with  a  given  set  of  noises  driving 
both  the  state  and  observation  vectors,  has  a  great  deal  of  descriptive  power 
as  well  as  mathematical  elegance. 

The  second  point  to  note  is  that  for  this  particular  type  of  model  an  exact 
state  density  solution  exists.  It  exists  in  a  number  of  forms,  one  of  which 
is  the  Kushner  equation  and  the  second  is  the  representation  theorem.  Since 
this  solution  exists,  one  can  in  theory  solve  for  an  optimal  tracker  under 
the  conditions  of  the  model  and  this  optimal  tracker  is  completely  mathematically 
wel 1  d  ef i ned . 

All  necessary  bounding  properties  of  a  given  tracker  can  be  extracted  from 
looking  at  the  properties  of  the  solution.  The  solution  is,  however,  computa¬ 
tionally  expensive,  and  statistical  as  well.  Monte  Carlo  runs  may  be  required 
to  find  the  average  properties  of  the  solution  and  Monte  Carlo  runs  of  the 
computationally  intensive  solution  can  be  prohibitively  expensive. 

By  studying  the  exact  state  density  solution  one  can  get  insight  into  the  proper 
ways  in  which  one  can  design  a  tracker.  As  an  example  of  this,  an  intelligent 
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tracker  design  recently  proposed  by  Hughes  was  motivated  from  a  study  of  the 
exact  solution  generated  by  the  Kushner  equation. 

Since  solving  the  Kushner  equation  in  real  time  is  impractical  with  the  current 
state  of  technology,  and  even  solving  it  off-line  as  a  method  of  analysis  can 
be  impractical  for  a  tracker  containing  a  large  number  of  states,  the  use  of 
analytic  techniques  to  bound  and  approximate  the  properties  of  the  solution  has 
been  examined.  Among  these  are:  the  Bobrovsky- Zakai  method  which  is  basically 
an  extension  of  previous  work  on  providing  a  Cramer-Rao  type  bound  on  the 
performance  of  a  system;  the  technique  of  cone  bounded  non-linearities  which  can 
bound  tracker  performance  relating  to  the  degree  of  non-linearity  of  the  system; 
extended  Kalman  techniques  which  can  provide  a  very  precise  indication  of  the 
accuracy  one  would  obtain  for  small  variations  in  target  state;  and 
methods  of  analyzing  non-linear  stability  in  terms  of  calculating  the  time 
for  a  given  tracker  to  completely  lose  track  of  the  target. 

At  this  point,  it  is  seen  that  several  extensions  and  improvements  to  these 
given  techniques  would  be  very  desirable.  Tightening  the  existing  bounds, 
experimenting  with  various  Liaponov  functions  for  computing  stability,  and 
extending  the  Bobrovsky-Zakai  bound  to  handle  Poisson  inputs  are  three  areas 
that  could  benefit  from  more  study. 

Once  one  has  the  analytic  techniques,  in  principle  one  can  es  ‘mate 
how  close  to  optimal  any  given  sub-optimal  design  will  fall,  and  one  can 
derive  the  functional  relationships  between  the  system  parameters  which  would 
be  of  particular  utility  in  the  specification  of  trackers. 

These  bounds  and  approximations  that  could  be  useful  for  estimating  and  specify¬ 
ing  the  performance  of  the  tracker  have  been  verified  directly  by  simulation. 

In  the  process  several  solutions  to  Kushner' s  equation  have  been  generated. 
Applications  have  been  explored,  specific  cases  of  tracking  N-state  random 
target  motion  have  been  studied  with  respect  to  jitter  properties  and  a  pro¬ 
cedure  by  which  one  can  describe  the  general  and  formal  procedure  for  defining 
and  designing  a  tracker  has  been  initiated. 
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Appendix  A  -  Application  of  MACSYMA 

MACSYMA  is  a  powerful  system  of  computer  programs,  developed  at  MIT, 
for  performing  computer  algebra.  (An  interesting  article  on  computer  algebra 
appears  in  the  December,  1981  Scientific  American.)  MACSYMA  has  the  capability 
of  symbolically  performing  algebra,  expansion  and  distribution  of  terms, 
combining  like  terms,  factoring  and  other  types  of  simplification,  calculus 
including  differentiation  and  integration,  solution  of  differential  equations, 
matrix  calculus  including  multiplication  of  matrices  and  finding  inverses, 
and  more. 

MACSYMA  has  been  used  in  the  past  to  perform  calculations  that  either  have 
or  would  have  required  months  or  years  to  perform  by  hand.  For  example,  com¬ 
puter  scientists,  in  about  20  hours  of  computer  time,  have  obtained  formulas 
which  reproduce  (and  correct)  formulas  first  obtained  by  the  18th  century 
astronomer,  Charles  Delaunay,  which  he  published  in  two  volumes  after  twenty 
years  of  derivation,  checking,  and  rechecking.  Surprisingly,  only  three  small 
errors  in  Delaunay's  work  were  found.  Further  work  has  modeled  effects  not 
treated  in  Delaunay's  work. 

This  points  out  the  main  use  of  MACSYMA  for  either  greatly  reducing  the 
time  required  to  derive  and  check  a  difficult  result,  or  in  obtaining  results 
which  one  could  not  obtain  at  all  by  hand. 

This  work  was  done  with  the  aid  of  MACSYMA,  a  large  symbolic  manipulation 
program  developed  at  the  MIT  Laboratory  for  Computer  Science  and  Support 
by  the  National  Aeronautics  and  Space  Administration  under  Grant  NSG  1323, 
by  the  Office  of  Naval  Research  under  Grant  N00014-77-C-0641,  by  the  U.S. 
Department  of  Energy  under  Grant  ET-78-C-0”  -4687,  and  by  the  U.S.  Air  Force 
under  Grant  F49620-79-C-020.  Special  thanks  are  due  Leo  P.  Harten  for 
supplying  valuable  information  on  use  of  the  MACSYMA  system,  and  for  helping 
us  overcome  problems  we  have  encountered  In  using  It. 


Our  Use  of  MACSYMA 


Because  we  were  continually  running  Into  problems  involving  lengthy, 
tedious  algebra,  we  decided  to  use  MACSYMA.  Our  initial  use  of  MACSYMA  was 
to  perform  algebra  involved  in  calculating  bounds  formulas  for 
translational  and  rotational  target  models.  Although  rigorous  formulas  were 
obtained,  the  bounds  proved  to  be  too  small  in  practice  to  be  of  practical 
use. 

Consequently,  an  approach  based  upon  linearizing  the  measurement  model 
was  undertaken.  This  led  to  matrix  Rlccati  equations  for  the  various  models: 
A  2x2  matrix  Rlccati  equation  for  the  one-dimensional  dynamic  model  with 
white  noise  driving  force;  a  3x3  matrix  Rlccati  equation  for  the  same  model 
with  Poisson  type  force;  a  4x4  matrix  Riccati  equation  and  a  6x6  matrix 
Riccati  equation,  which  generalize  the  previous  two  models  to  permit  y-motion 
as  well  as  x-motion. 

Even  the  simplest  of  these  would  be  extremely  difficult  to  do  by  hand, 
involving  tedious  checking  and  rechecking.  However,  with  the  use  of  MACSYMA, 
closed  form  expressions  have  been  obtained  for  all  of  these  problems  and 
verified,  an  extremely  difficult  task  to  do  by  hand,  since  an  explosion  of 
terms  occurs  before  final  cancellation. 

MACSYMA  has  also  been  employed  to  derive  the  form  of  the  differential 
equations  which  must  be  used  to  obtain  the  Rhodes-Gilman  cone-bounded  non¬ 
linearity  bounds. 

Computer  printout  of  some  of  this  work  is  given  here.  Five  batch  programs 
were  used  and  1  Is  given  here.  In  a  batch  program  the  instructions  are 
read  from  a  previously  prepared  file  and  executed  one  after  the  other.  The 
Instructions  appear  on  the  C  lines  in  lower  case,  while  MACSYMA’s  responses 
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appear  on  the  0  lines.  Numerous  explanatory  conments  are  given  and  have  the 
form: 

I*  This  Is  a  connent  */. 

Batch  program  A  shows  the  derivation,  solution,  and  verification  of  the  2x2 
matrix  Riccati  equation  associated  with  the  1-D  dynamical  translation  white- 
noise  force  problem.  The  matrix  Riccati  equation  appears  in  lines  (C62)  and 
(D62) .  The  solutions  are  given  by  (055),  (D56),  (057),  and  (D58)  where 
D(t)  is  defined  by  (D48).  Derivation  is  given  by  the  preceding  steps, 
while  verification  is  given  in  subsequent  steps.  The  matrices  A,  M,  and  N 
used  in  defining  the  Riccati  equation  are  defined  at  the  beginning  of  the 
printout.  Note  that  the  form  of  the  Riccati  equation  assumes  that  the  time  is 
normalized  in  a  suitable  way  (with  no  loss  of  generality).  Note  the  symmetry  of 
the  matrix  elements  as  required.  A  rough  computer  plot  of  the  matrix  elements 
is  also  given. 
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(02)  0atch(ric22t>,>l4«k,«ilcOK)l 


(Cl J)  /*  fhlni  diffiriitul  *ii.  far  l*! 


i«li<ilf(i,(l  *  b.i; 
lint*  712  nut. 


(Ill) 


(  4 

I  --  (Xll(»>) 
(  41 
C 

t  d 

t  —  (Midi) 
(  41 
t 

I  4 

t  --  (YIK1I) 
l  41 
I 

(  4 

c  --  mi  (in 

(  41 


4 

—  1X12(1)1 
4T 

4 

—  (122(1)1 
41 

4 

—  1112(1)1 

41 

4 

--  1122(D) 

41 


(  121(11 
t 

t  f  121(11 
[ 

•  (  XIK1I 

I  - 

1  F 
t 

t  -  111(1) 


ICI4)  M»U* 
lint*  4  Mic. 

(CIS)  for  id  thru  4  4o  (for  jtl  thru  2  4o  •4»iciM(yirt<z»4,M,)>»p»ft(zn,2f  I, 
lint*  142  mic. 

(CIA)  /»  iMcily  wikioiM  t*  b#  tilrid  for  •/ 


until 

llN>  4  MIC. 

(C17)  fir  ill  thru  4  do  (fir  jil  thru  2  do  uncoitlzli,  j),ui>>4 
I  in*  112  MIC. 

(CIO  /•  lit  iiitiol  coiditioni  •/ 


xj»itiu*<t(0,t.x>aiiittri«(2f  2,1,1,11; 
liM*  II  MIC. 


[  XII 40) 

XI 2 ( 0  >  1 

(  0 

1  1 

(III) 

( 

1  • 

( 

J 

t  X21I0I 

122(0)  1 

(  0 

*  1 

(«♦)  yiiititu»it(«,t,y)ai4iit<2); 

IlM*  11  MIC. 

[  11110) 

112(0)  J 

(  1 

•  1 

(Ilf) 

( 

J  • 

( 

2 

(  121(0) 

122(0)  1 

(  0 

1  1 

(Ml)  Ur  j„  tfcr*  j  4#  (for  Jil  thru  2  4o  itriliodll, 


122(1) 
122(11 
XI 241 > 

F 

■  112(1) 


,iqi)  II 
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(C21)  for  i 1 1  thru  2  do  tt or  jit  thru  2  da  aUilue(y[ilj},t:>0,pii'tl/i>it(2,il Jill! 

1 1 «•=  83  nsec. 

Tine*  342?  nsec. 

(C23)  batch! ric22c , >,dsk,ui lcox) t 

(C24)  /*  Solve  lifferential  Eqs.  for  X  and  Y  */ 

illlldlUl«t(ltS,Ul)l 
fine3  SI1S  nsec. 

IC25)  »nniub»t<<lv»r‘2«Zii»<lvar*2-Zi),lvar'4*1,ao»  >« 
line*  123  nsec. 

IC242  for  til  thru  I  do  aMUiJtpart<ano(i)>l)*ilt<parlfanoCi)'2,12llv»r,t)t 
Tine*  3432  nsec. 

IC22)  for  ill  thru  8  do  »n«1  £i  liratexpandtntiubnt  <  t/»qrt  (2  2  ,«qrt  <2)*t/2,r»alperUan«U  1  2  2  2  )  I 
Tine*  24221  nsec. 

(C28)  for  ill  thru  I  da  pri*tt»Mt(i32t 
T  I 

Y22IT)  *  C0S( - >  COSHt - > 

SQRT<2)  SORT  <  2 ) 

IT  T  T 

COSt . 1  SIHHt . )  COSHt . —  1  SINt . 2 

SORT  f  2  >  SORT (2)  SORT (2)  SBRTI2) 

Y2I(T)  * . . . - . 

SORT (2)  SORT (2) 

T  T  T  T 

COSHt . >  SINt . 2  cost . 2  SIHHt . 2 

SORT (2)  SORT  <  2  >  SORT  <22  S0RTI22 

TI2IT)  * - - 

SORT  121  SUfc  T  <  2 1 

I  T 

TM(I)  *  COSt . 2  COSHt - 2 

SOOT <2 1  SORT  12) 

T  T  T  T 

F  COSt . )  SIHHt . I  F  COSHt - 2  SINt . 2 

SBRK2)  SORT  1 2  2  SORK2I  S08K2> 

X22IT 2  * . —  *  — . . . 

SORT  <  2  2  SORT (2) 

T  T 

X2KT)  *  -  F  SINt - 1  SIHHt - > 

SBRTI2I  SORT (2) 

T  T 

X)2tT )  *  F  SINt - 2  SIHHT - 2 

SORT  122  SORT  1 2) 
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IT  T  T 

f  COSI . I  SIMM l - 1  f  COSHI - )  SiN< - > 

SURTI2)  SOS  T  ( 2  >  SORTI2)  SUITI2) 


S0RTI2) 


SORT ( 2  > 


XIHT)  * . 

Tint*  (42  nttc. 

(C2f)  /•  Obtain  X  and  T  natricts  »/ 

xiiubttdaat  1(91,  anti  C4)lantU2IlantK83]lx>l 
lint*  SO  nine. 


(CTO)  yifub«t((an<1(1),anft[2),aMH3]|aan1C4]]>y)t 
Tint*  51  nttc. 

(C1I)  print!"!  •  *,«)» 

T  I  T  T 

f  COS! . )  SINN! . )  f  COSH! . )  SIN! . > 

SORT! 2)  SORT  12)  SORTI2I  SORTI2) 


X  * 


SORT  1 2) 


SORT! 2) 


T  T 

f  SIN! - )  SIMM . ) 

SORT (2)  S0RTI2) 


I 


f  COS! - 


T  T 

-  r  SIN! . )  SINH! . ) 

SORT  12)  SORT ( 2 ) 


SORT  12) 


Tint*  4IS  nttc. 

IC32)  print!"!  *  ",y)l 


T  * 


T  T 

COS! . )  COSH! . I 

SORT ( 2 )  SORT! 2) 

T  I  T  T 

COS! . )  SINH! . I  COSHI . I  SIN! . ) 

SORT (2)  SORT! 21  SORT  <  2  >  S0RTI2) 


I  T 

I  SINH! - )  f  COSHI - )  SIM C  — . ) 

SORT (2)  SORT! 2)  SORII2) 

- - - 4 - - - - 

SORT (2)  SORTI2) 


T  T  T  T 

COSH! . I  SIN! . )  COSI . )  SINN! . — ) 

SORT  12)  SORT! 21  SORTI2)  SORT! 2) 


SORT  12) 


SORT  12) 


SORTI2) 


$001(2) 


Tint*  ITS  nttc. 
lint*  41343  nttc. 


I  T 

COS! . )  COSH! . ) 

SORT (2)  SORT  12) 


IC34)  »atcMric22d,>,dik,»ilcax)* 
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(Cli)  /»  Obtain  inverse  of  I  >/ 

Oetautitrual 
Tine*  )  nsec . 

ICJ4)  dtallnxtps ifnlst* 

I  let*  }  estc. 

(C17I  yinviy‘*-lt 
line*  27 J  nsec. 

(CJI)  /•  Obtain  Otltrniaait  Dtnuninator  •/ 

•Jtl;4enon(yiavll 
lint*  1  nsec. 

(CIO I  /*  Itfite  substitutioss  •»/ 

subtisiiiT>H/i'i-t<  2)  >*(t>pU/t<|rt(2>>-txp(-l/t')rt  121)1/2; 

Tint*  741  nstc. 

T  T 


SORT  ( 2 )  SQRTI2) 

T  Zl  -  XI 

( 049)  SIMM! - - - — . . 

SORT  12)  2 

(CTO)  sut2:coth(t/tqrl(2))*(txp(t/i>]rt(2>)*fxp(-t/sorl<2>>>/2; 
lint*  104  nstc. 

T  T 


S0RTI2)  SORII2I 

T  «  ♦  XI 

TR401  COSHI— . . . - . . 

SORT  12)  2 

IC4I )  sutlli(sinlt/si)rt(2)))‘2; 

Tint*  124  nstc. 

2  T 

(HI)  SIN  I - ) 

SORT (2) 

(C42)  sub]ril-(c«sll/t4i't(2)l)  2; 
lint*  121  nsec. 

2  I 

(M2)  1  -  COS  I . ) 

SIRII2) 
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(C4J)  iub4lla *(2»t/Mrt(2)>*(c*»h<t/l<|rt(2>M*inb<t/iqrt(2)l>"2; 
Tim*  221  mc< 

2  T 


S80K2)  1  12 

I Mil  U  ■  (SINN! - )  4  C8SHI - II 

(OKI 121  8001(2) 

IC44I  iub2iIa‘<-2*t/aqrt(2>>«(ca>Mt/iqrt<2>)-fink(t/iqrt<2))>‘2; 
lint*  222  M*c. 

2  1 


8001(21  T  T  2 

(8441  U  ■  ICKIK - 1  -  tINHt - II 

8801(21  88*1(21 

8C4SI  »ub4li(»iab(t/«qrt(2l))*2; 
liaa*  H  kmc. 

2  1 

11491  SIMM  I - 1 

8801(21 

IC44I  tabbri (cum t/aqrt(2l ll'2-l; 

IlM*  fl  MM. 

2  1 

<8441  COSH  < - 1  -  I 

8801(21 

IC47I  /•  Hat*  »*•»•  aubatiluliaai  la  Sataraiaaat  •/ 

•)aairatiabat<luMrlUib41lrat*xpaa4<«*bat(«ab4>iukal<tub9lrtta!<8ao0<'‘*l«iMt<*iiblrlaabl],r*taxpka4(iub*l((t*0lftub2]l4at)ll)))))8 
1 1«»*  1989  a*M. 

4 C48 1  4*414(11  •  4ta; 

Iim*  9  a*M. 

2  1  2  1 

<8481  8(11  •  COW  < . I  ♦  C88  I - 1 

8801(2)  8081(2) 

(C48)  /•  Obtaia  Huaaratar  af  1  lavtraa  */ 

ylavauaiaualyiavlt 
TIm*  9  mm. 

IC90I  /•  8btala  Huaaratar  of  *  •  I  .  (1  iawraal  */ 
vauaia.ytaaauat 
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(C5I)  /•  Sinplify  Vmmtrator  usmi  -;•«*  substitutions  uu>l  to  sinpllfy  III)  »/ 

« 

(CSI)  for  ul  thru  2  da  (far  j:t  thru  2  da  ynunti , j  1  < rila.ipandlntsubsK tub3r, tubJI, rattxpoadf subst (Csubl ,sub2)fvnun( i , jl > ) ( > ) I 
lint*  1*31?  tut. 

tC52»  far  ul  thru  2  do  (far  jil  thru  2  do  viuw(t  ,.j  )i  rntsubsttsubor, sub*) ,  ratexpnndl  subst  < subt, subst (subS.vnunti,  jll  1 1 > It 
'I tat*  KSS  nstc. 

(CS3I  /«  Tha  solution  to  tt>«  2*2  lUtrin  Kiccali  Equation  is  nun  qinan  by  */ 

yinatrix((vltttllyt2(U),lv2m>.v22(tl  II* 
linn-  ID  him. 


(C5d)  printlaV  •  *,«lt 
(  Vtl(f)  V12I1I  ) 
V  *  l  3 

(  V2IIII  V22ITI  1 
I mt2  $1  ante. 


(CSS)  /•  uhti-t  Villi), Vl2lt),V21(l),V22(l)  ora  jivM  by  */ 

tqHlvll(t)*VMlnC1,l]/d(t); 
lint*  15  nsac. 

IT  IT 

2  f  CQSH( . >  SIMH( . )  -  If  COtl . I  SIN( . I 

SOR 1 1 2 >  suit  112)  sum  (2)  S0RH2) 

(D55>  VI I  (II  *  . . - . 


$1)11121  KT1 


(C5a)  tql 2:«I2( t )-vnun( I , 2)/d( t I ; 
lint2  15  me. 


<  03o) 


V12II) 


2  T  2  T 

r  cosh  ( . >  -  r  COS  I-- . I 

SH)I(?)  SQITI2I 


1(1) 


ICS?)  »i|2l  1*21  <  t  )“*nun(2,  IJ/dlU ; 
lint2  15  attc. 


(05?) 


V2MI)  2 


2  T  2  1 

f  COSH  ( . I  -  f  COS  ( . I 

S0RK2I  SOItt  (2) 


1(1) 
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(CM)  eq22:v22(tlsvaun£2,23/d(t ),' 
ItM*  IS  HSSC. 


USD 


V22( 1 )  * 


T  T  t 

2  F  COSH ( - )  SINH( - )  ♦  2  f  C0$< . 

SORT  <  2 1  SURK2I  SQRK2) 


SUR 1(2)  D ( T ) 


(C5»)  /»  uhere  Kt)  it  defined  above  */ 

/*  Rote  that  I)  ie  a  synnetric  natrix  at  required  */ 


line*  32 set  atec . 

(CatO )  batch(ric22e,>,4sk,uilcox)t 

(Ctl)  /*  Verify  that  V  satisfies  the  2«2  Riccati  equatioe  •/ 

/•  First  substitute  VI2(t>  for  V2l(t)  in  U  natrix  •/ 

printCV  x  •,visubst(vl2(t),v2l(t),v))S 
£  VII  IT)  V12IT)  3 
V  «  C  ] 

(  V 1 2 ( T )  V22IT)  J 
line*  58  nsec. 


(Ct2)  nceqidif f ( v, t )«a.v*v .transpose! i)*n-v.n.v; 
line*  Ml  nsec. 

[  2 

t  d  d  }  t  VII  (T)  „ 

t  —  (VIKTI)  —  (V12III)  J  £  2  V12IT) .  V22(I) 

£  41  4T  3  £  F 

(It2)  £  3  »  l 

£  4  4  3  t 

E  —  (VI2(D)  —  (V22(T1)  3  £  VI Ml)  V12IT) 

£  4T  41  3  £  V22(l>  -  F 

E  F 


(Ctl)  ricif 3t 
line*  4  nsec. 

(C44)  for  til  thru  2  4o  (for  jii  thru  2  4o  ric icens(part(riceq,l , i, j)*part(riceq,2, i , j), 
line*  22  nsec. 


T 

)  Slid . — ) 

SORT (21 


VI I ( T)  VI2II) 
F 

2 

VI2  (T) 

F 

ric))f 
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(C6j)  for  ill  thru  3  do  pnntlracCi  lit 
2 

J  1112  III 

-  (V22(T>)  *  F . 

41  f 

•J  VI I II I  VI2ITI 

—  IVI2IIII  *  V22UI  -  — . 

41  f 

2 

J  VII  III 

-  (Vllllll  *  2  VI2ITI . — 

.Jl  f 

Tim*  200  met. 

ICOO)  /*  Oifforoiitiato  the  Solution  */ 

<>ol*iCeo22l*4l2i*4nil 
line-  5  nsec. 

(Co2)  JsolniJiff (soln,t)» 
lint-  1042  mat. 

IC08  >  for  111  thru  3  Jo  print < Jsolnl i 3 )» 

2  1  2  1  2  1  2  1 

2  f  SUN  I . I  2  F  SIN  I . I  2  F  COSN  I . >  2  F  COS  I . ) 

SOK I (21  SORT  1 2)  SORT  12)  S0RI12) 

-  -  - . . * -  t  - 

j  S0RTI2I  SURII2I  $aRT<2>  SORT  12) 

-  IV22II))  =  . - . . — 

Jl  SORT  1 2)  0(1) 

II  114 

(2  f  cusHi— • — “r*si«hi - )  *  2  r  cost . -»  sui . >>  <--  (»n>)> 

S0RK2)  S0RK2)  SQRII2)  SUN  1(2)  Jl 


SORT ( 2 )  0  (I) 

I  I  I  •  I 

2  f  COSHI . I  SINIK . )  2  F  cost - )  Slid . ) 

SURI  (2)  SOM  (2)  SUR  1(2)  SURTI2)  2  1  2  14 

.  * .  (F  COSH  1 - )  -  f  COS  1 . ))  1“  10(1))) 

J  SORT (2 1  SBRTI2)  SORT (2)  SORT (2)  41 

-  (Vl2(l>)  *  . - . . - . - . 

Jl  0(1)  2 

0  (I) 
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2  1  2  1  2  1  2  1 

2  F  SIHH  ( . I  2  F  SI*  ( . )  2  F  COSH  ( - 1  2  F  COS  I— . ) 

SORT (2)  SORT (2)  SORT (2)  SORT(2> 

-  ( - *  -  - - 

<1  SORT!  2)  SORTI2)  SORT  (2)  SGRTI2) 

—  I VI I  IT ) )  t  - - 

dT  SORT  12)  BIT) 

T  T  T  T  d 

(2  F  COSHI . I  SINH1 . I  -  2  F  COSI . — )  SI*t . ))  I--  (PIT))) 

SORT < 2 >  SORT (2)  SORT  1 2 )  SQR T 1 2 >  dT 


Tin**  1141  node. 

IC69I  /*  Verification  pr oceeds  bv  subtracting  above  set  of  equations  fron  these  »/ 

/*  and  noting  appropriate  substitutions  and  sinplifications  •/ 

for  ill  thru  3  do  rlcti]isub*t(ssla,dsoln[11-rlcCi)>t 
Tine*  112  nsec. 


2 

S0RTI2)  0  (T) 


(C20I  /*  Obtain  derivative  of  deterninant  «/ 

dddlidif f (dnq.l); 

Tine*  181  nsec. 

IT  T  T 

2  COSHI . )  SIHHI - >  2  COSI . )  SIHI . ) 

d  SORT  12)  SORT! 2)  SORT  1 2)  SOR  T  <  2  > 

1020)  -  I  BIT  > )  *  . - . - . 

dT  SORT  12)  SORT! 2) 


IC7I)  for  til  thru  3  do  rlctl]iratsubsllsubdr,subdlrratsubsti(ub3r,sub3l(rate>p>adlsubst(dnq,rats>np( <d< t) )'2*subst(dddt,rict i Jl ) ) I ) 
! 

It 

Tine*  7531  nsec. 

(C72)  for  ill  thru  3  do  printlrictillt 
0*0 
0  <  0 
0  *  0 

Tine*  *0  nsec. 

IC73)  /*  This  conpletes  the  verification  */ 


Tine*  13033  nsec. 
Tine*  fi8}9  nsec. 
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Appendix  B  -  Representation  of  Rectangular  Solids  by  Vertex  Projections 
The  target  model  is  assumed  to  be  a  rectangular  solid.  For  the 
tracking  problem,  it  is  necessary  to  determine  the  silhouette  of  the 
target  projected  onto  the  image  plane.  Figure  1  defines  the  notation 
for  the  target  model . 


The  x’,  y' ,  and  z‘  axes  are  fixed  relative  to  the  target  with  the 
original  coinciding  with  the  center  of  mass.  The  front  of  the  tank  points 
in  the  direction  of  the  x'  axis,  while  z'  is  vertical  and  «'  points  off  to 
the  left.  It  is  assumed  that  at  t  =  0,  the  x',  y' ,  and  z'  axes  coincide 
with  earth-fixed  x,  y,  and  z  axes,  but  that  in  general  they  will  be  trans¬ 
lated  and  rotated  with  respect  to  each  other.  To  simplify  the  analysis,  it 
is  assumed  that  z  =  0  at  all  times.  That  is,  the  tank  moves  only  on  level 
ground.  However,  the  position  of  the  tracker  with  respect  to  the  target  is 
arbitrary.  In  general,  the  tank  makes  an  angle  9  with  respect  to  its  initial 
direction  as  shown  in  Figure  2. 

x' 


Figure  2.  Orientation  of  Target  as  Seen  From  Above 


Thus,  the  fixed  and  moving  coordinate  axes  are  related  by 
x  =  xc  +  x'  cos  e  -  y'  sin  0 
y  =  yc  +  x'  sin  e  +  y*  cos  0 
z  =  z' 

where  xc  and  yc  are  the  x  and  y  coordinates  of  the  centroid. 

While  the  position  of  the  target  is  assumed  to  be  given  by 
r  =  xi  +  yj  +  zk, 

where  i,  j,  and  k  are  an  orthogonal  set  of  unit  vectors,  the  position  of 
the  tracker  is  assumed  to  be  given  by 
R  =  Xi  +  Yj  +  Zk  . 

The  unit  vector  n  defining  the  direction  of  the  target  as  seen  by  the 
tracker  is  given  by 


We  specialize  further  to  the  situation  where  the  tracker  and  target  are 
far  from  each  other  so  that  n  is  constant  for  all  practical  purposes  as  the 
target  goes  through  its  motion.  In  effect  then,  only  the  angular  orientation 
affects  the  projected  silhouette  on  the  image  plane. 

In  addition  to  n,  two  orthogonal  unit  vectors  p  and  q  are  defined  by 
p  x  n  =  q  . 

p  and  q  define  the  direction  of  the  "horizontal"  and  "vertical"  axes  for  the 
tracker  image  plane. 

The  projection  for  each  of  the  eight  vertices  onto  the  image  plane  is 
required.  For  any  arbitrary  vector  r,  this  projection  Pr  is  given  by 

Pr  =  (p«r)P  +  (q-r)q 

A  A 

=  PP  +  qq  . 
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The  unit  vectors  n,  P,  and  q  are  related  to  the  unit  vectors  i,  j, 

and  k  by 

n  =  j  cos  tp  -  k  sin  \p 
P  =  i 

q  =  P  x  n  =  k  cos  ip  +  j  sin  ip 
or 

i  =  P 

j  =  n  cos  ^  +  q  sin  <p 

k  =  q  cos  i|»  -  n  sin  ip 

where  tp  is  the  depression  angle  of  the  target  as  seen  from  the  tracker. 
ip  tan  y 

Thus,  the  "horizontal"  and  "vertical"  components  in  the  image  plane,  p  and  q, 
are  given  by 

p  =  p-r  =  i  -r  =  x 
q  =  q-r  =  y  sin  ip  +  z  cos  i|> 

Coordinates  of  x' ,  y'  and  z1  are  as  follows: 


Point  § 

x' 

y' 

z' 

1 

-a/2 

-b/2 

-c/2 

2 

a/2 

-b/2 

-c/2 

3 

a/2 

-b/2 

C/2 

4 

-a/2 

-b/2 

c/2 

5 

-a/2 

b/2 

-c/2 

6 

a/2 

b/2 

-c/2 

7 

a/2 

b/2 

c/2 

8  • 

-a/2 

b/2 

c/2 
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From  these  coordinates,  the  projections  onto  the  image  plane  are  as 


follows: 


Pi 

= 

xc  - 

H  a 

cos 

e 

+  % 

b 

sin 

0 

p2 

= 

xc  +  %  a 

cos 

e 

+  % 

b 

sin 

0 

P3 

= 

p2 

P4 

= 

pl 

P5 

= 

xc  ' 

H  a 

cos 

0 

-  >5 

b 

sin 

0 

p6 

= 

xc  " 

■  H  a 

cos 

e 

-  H 

b 

sin 

0 

p7 

= 

p6 

p8 

= 

P5 

*>1 

= 

sin 

*  *c 

-  *5 

a 

sin 

6 

-  *5 

b 

cos 

0 

-He 

cos 

* 

q2 

= 

sin 

*  yc 

+ 

a 

sin 

0 

-  % 

b 

cos 

0 

-  H  C 

cos 

<3 

= 

sin 

*  yc 

+  % 

a 

sin 

0 

-  H 

b 

cos 

0 

+  H  c 

cos 

* 

<*4 

= 

sin 

*  yc 

-  H 

a 

sin 

e 

-  H 

b 

cos 

0 

+  H  c 

cos 

* 

q5 

= 

sin 

*yc 

-  >5 

a 

sin 

0 

+  % 

b 

cos 

0 

“  H  c 

cos 

* 

q6 

= 

sin 

*  yc 

+  *5 

a 

sin 

0 

+  *5 

b 

cos 

0 

-  H  c 

cos 

q7 

= 

sin 

*  yc 

+  ’s 

a 

sin 

0 

+  >s 

b 

cos 

0 

+  H  c 

cos 

>!' 

q8 

= 

sin 

*  yc 

-  4 

a 

sin 

e 

+  h 

b 

cos 

0 

+  H  c 

cos 

<l> 

Now  it  may  be  shown  that  the  quantity 

(q  -  qjMPj  -  Pf)  -  (p  -  PfMqj  -  q,.) 

is  zero  on  the  straight  line  between  the  two  points  i  and  j,  positive  in 
the  region  to  the  left  of  the  line  in  going  from  i  to  j,  and  negative  in 
the  region  to  the  right  of  the  line  in  going  from  i  to  j.  (It  is  also 
easily  shown  that  this  quantity  is  invariant  to  coordinate  translations  and 
rotations.)  Hence,  the  target  silhouette  can  be  defined  in  terms  of  the 
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following  function 


=  [}..  (p.q)  =  U((q  -  qi ) ( Pj  -  pi )  -  (p  -  Pi ) ( qj  -  q^) 

where  U(x)  is  the  unit  step  function 
i 1,  x  >  0 

U(x)  = 

f  0,  x  <  0 

A  different  projection  (p.q)  is  obtained  depending  upon  which  corner 
number  k  is  closest  to  the  tracker.  One  finds 

u(3)  =  U12  U26  U67  U78  U84  U41 

U(4)  (P-P)  =  U5i  U12  U23  U37  U78  U85 

U(?)  =  U26  U65  U58  U84  U43  U32 

U(8)  (p.q)  =  Ug5  U51  Uu  U43  U3?  U?6 

where  the  dependence  upon  p  and  q  has  been  suppressed.  (These  results  have 
been  derived  with  the  aid  of  a  three-dimensional  solid  rectangular  model 
whose  corners  are  numbered.) 

These  four  quantities  can  be  thrown  into  a  single  formula  for  the  pro¬ 
jection  of  the  target  on  the  image  plane,  valid  for  any  orientation  so  long 
as 

0  <  ip  <  90°, 

J(p.q)  -  U(-o)  U(  8)  U(3)  (P.q)  +  U(-a)  U(-B)  U(7)  (p.q) 

+  l)( a)  U(-6)  U(8)  (p.q)  +  U(a)  U( 6)  U(4)  (p.q) 
i  J(p,q;  Xc.  Yc,  e,.  v) 

where  a  and  0  are  the  direction  cosines  defined  by 

a  =  n-i' 

6  *  n-j' 
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Since 


where 


+  (yc  +  R  cos  ty)  j  -  R  sin  ^  k 


and 

A  A  A 

i '  =  1  cos  e  +  j  sin  e 
j'  =  j  cos  e  -  i  sin  e 

4  4  a 

a*  =  (rt  -  R)  ■  i '  =  xc  cos  e  +  (yc  +  R  cos  i|»)  sin  0 

6'  =  (?c  -  R)  •  j'  =  (yc  +  R  cos  tl>)  cos  0  -  xc  sin  0 

Since  a1  and  a'  have  the  same  sign  as  a  and  @  ,  respectively,  one  can  also 

» 

write  U(p,q;  xc,  yc<  0,  *) 

*  U  (-a*)  U(B')  U(3)  (p,q)  +  U(-a')  U(-B’)  U(3)  (p.q) 

+  U(a')  U(-e')  U(8)  (p,q)  +  U(a')  U(b')  U(4)  (p,q)  . 


This  result  Is  valid  for 
0  <  *  <  90°. 

For  the  case  where  <(/  =  90 
U(p,q;  xc,  yc. 


°,  i.e.,  looking  straight  down  on  the  targets 
0,  90°)  =  U43  U37  U?8  U84  . 
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APPENDIX  C  WHITE  GAUSSIAN  FORCING  FUNCTION 


This  appendix  fills  in  the  details  of  the  calculation  for  the  1-D 
translational-motion-only  white-noise  force  model.  One  begins  by  converting 


the  basic  dynamical  equation 


mx  =  F$  +  F  ( t)  -  Kx 


into  a  first-order  system  in  terms  of  state  variables  x^  and  x2  defined 

relative  to  the  average  velocity  F  /k  produced  by  the  steady  force  F 

s  s 

as  follows: 

x1  =  x  -  (Fs/K)t 
x2  =  X1  Fs/K 


X1  =  x2 


x2  =  -  Dx2  +  ar(t) 


where  the  damping  constant  D  is  defined  by 

0  =  k/m 

and  the  random  acceleration  a  (t)  is  defined  by 

r 

ar(t)  =  Ff  (t)/m 


where  Fr ( t )  is  the  random  white-noise  force.  Combining  the  state  variables 


x.  and  x?  into  the  state  vector  X  defined  by 


the  stochastic  dynamical  equation  for  the  system  is  given  by 


dx  =  A~x  dt  +  B  dw, 

ss  ~  I 


where  the  dynamical  matrix  A  is  defined  by 


!o  l\ 

\°  ’D/ 


- r  i 


B  Is  defined  by 


B 


and  ^wr  is  defined  by 


■C) 


dwr  -  a^(t)  dt 


Since  ar ( t )  is  a  white-noise  process,  vr(t)  is  a  Brownian  motion  process 
satisfying 

<(dwrdwr)>  =  Q  dt 


This  completes  the  statistical-dynamical  specification  of  the  state 
variables.  To  proceed  further,  it  is  necessary  to  specify  the  measurement 
process.  The  measurements  are  all  assumed  to  be  contained  in  an  n-component 
column  vector  g“  (x)  /g ^(x)\ 

92m\ 

?<x>  *  ' 


\gn(x) / 

X  L 

where  9^(x)  is  the  temperature  measured  by  the  ixn  pixel  (in  the  absence  of 
measurement  noise)  and  it  is  assumed  that  the  sensor  has  n  pixels  in  all. 
Defining  z  as  an  n-component  column  vector  whose  components  represent 
Brownian  motions,  i.e.dz  represents  a  white-noise  increment  added  in  a 
time  dt,  the  total  observed  increment  to  the  measurement  process  is 

9iven  b*  dy  =  “g  (7)  dt  +  dT 


It  is  assumed  that  the  noise  process  for  different  pixels  are  uncorrelated, 

so  that  the  noise  spectral -density  matrix  R  is  given  by 

R  =  R  I 

as  ps 

where  |  is  the  n  x  n  idenity  matrix,  and  R  is  the  spectral  density  of  any  pixel, 
1  .e. 

<(dZ1)2>  =  R  dt 
C-2 


for  any  pixel  i. 

The  quantity  g^x)  is  defined  by  g.{x)=e.(x)  aT  where  0.(x)  is 
the  fraction  of  the  pixel  area  which  is  illuminated  by  the  target's  image. 

If  the  i**1  pixel  is  completely  inside  the  target's  image  when  the  target  is 
at  x  then 

e^x)  =  1, 

while  if  the  pixel  is  completely  outside 

0.(x)  =  0 

Now  the  general  formalism  involved  linearizing  the  measurement  vector 
about  the  estimate  of  the  state  vector,  leads  to  a  Riccati  equation  of  the 
form 

V  =  AV  +  VAt  +  M-  VN  V 

ss  ss  ss  ss  a  a  a 

where  A  has  been  previously  defined,  M  is  defined  by 

M  =  QBBT 

SS 


and  N  is  defined  by 


’12 


'll 


=  n21  =  n22 

ag  (x)  T  . 

R 

ax’ 


=  o 

ag  (x) 

■aT 


Since  ag^xj/ax  =  0  if  the  pixel  i  is  not  at  a  vertical  boundary, 
and  is  equal  to  ±  &T/e  if  it  is,  and  since  there  are  2  w/6  such  pixels  it 
follows  that 


N 


11 


1 

R 
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consequently 


.  ■  n 

\0  0/ 


where 


=  2  W 


53  R 


7  ^  (SNR)2  B 


It  simplifies  the  solution  some  and  gives  some  insight  into  how  the 
relaxation  time  depends  upon  various  parameters  to  convert  to  a  dimensionless 
time  variable  equal  to 

(Q  q)1^4  times  the  regular  time  variable.  When  this  is  done,  the 
Riccati  equation  has  the  same  form  except  that  the  definitions  of  the  A,  n,  and  N 

~  ~  B£ 

matrices 


A 

8 


M 


■( 

■( 


N  = 


0 

0 

0 

1/F 


1  1 
D' 


i) 


/,/h l  u\ 
Vo  o  / 


where 


D’  =  D/(Q  q)1/4 


F1  =  ( Q/ q3) 1 /4 
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and 


Under  usual  circumstances  O'  is  small  compared  with  1,  and 


therefore  may  be  neglected. 

The  solution  of  the  Riccati  equation  proceeds  by  solving  the  linear 
system 

X  =  A  X  +  H  Y 

Y  =  N  X  -  AT  Y 
~  x  «  »  ~ 

subject  to  the  Initial  conditions 

X(o)  =  0 

l(o)  -  I 

It  may  then  be  shown  that  a  V  defined  by 

V  =  X  Y'1 

~  as 

satifies  the  given  Riccati  equation. 

A  MACSYMA  program  to  solve  this  Riccati  equation  is  given  in  Appendix  A. 
(See  batch  program  A,  which  contains  explanatory  comments  for  the  derivation  and 
verification.  The  solution  is  given  by  D55,  D56,  D57,  and  D58  together  with 
the  definition  048.) 

Since  as  t  +  * 

v„./TT, 

and 

°TRACK  =  ^1/4  /F- 

j 

=  (4  Q/q3)1/8 
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APPENDIX  D- GENERALIZED  POISSON  FORCING  FUNCTION 


In  this  appendix  it  is  shown  how  the  method  of  linearization  is  used  to 
approximate  tracker  RMS  error  when  the  1-D  rectangular  target  defined  in  section  5 
which  satisfies  the  1-D  dynamic  equation 


mx  =  F  -  Kx 


(D-l) 


where  x  is  the  target's  location,  k  is  a  damping  constant,  and  F  is  the 
random  force.  The  random  force,  F  ,  is  assumed  to  consist  of  a  d.c.  term  F$ 
and  a  zero-mean  generalized  Poisson  term.  That  is. 


F  =  F$  +  Fr(t) 


(D-2) 


where  F$  is  the  d.c.  term  and  (t )  is  the  generalized  poisson  term  which  under- 

O 

goes  random  step  changes  of  force  with  mean-square  amplitude  oj:  and  frequency 

v 


Written  as  a  first  order  system,  x.  =x,  s?  =  x,  x,  =  x 

J  L.  j 


x}  *  x2 


k2  "  as  -  0x2  +  *3 


where  a$  *  Fs/m 


(D-3) 


af  s  Fr/m 
D  «  K/m 

The  quantity  Xj  ■  ar  is  a  generalized  Poisson  acceleration  which  undergoes 
a  Brownian  motion-type  amplitude  such  that  its  spectral  amplitude,  Q,  defined 


<(ar(t))2>=  Qt 


satisfies 


V  v 


%  2  ,2 

Xr  oF^/m 


(D-4) 


(D-5) 


2  2  2 

where  a,  -  at  /m  (the  mean-square  amplitude  of  the  generalized  Poisson 
«r  rr 

acceleration). 

From  equation  (D-3),  the  target's  stochastic  differential  equation  is 

d“x  = 

where 


A 


B 

ts 

and  ar(t)  is  the  generalized  Poisson  acceleration  with  spectral  amplitude  Q 
given  in  equation  D-5. 

The  measurement  is 


A  X  dt  +  B  d  a„  (D-6) 

fit  55  r 


dy  =  U(x)  dt  +  dv 


(D-8 


where 


g(x)  is  the  noise  free  observation  and 

E  {d"v  d  7T}  =  R  dt  .  (D-9) 

» 

Linearization  of  (D-8)  results  in  the  following  Riccati  equation, 

V  =  A  V  +  V  AT  +  a  -  V  MT  R"1  M  V  (D-10) 

s  gs  ss  st  m  9  s  »  ss  t9  a 
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where 


/O  0  0\ 

a  =10  0  0  J  (D-1 1 ) 

\0  0  p/ 

p  -  v  4  /m 2  (D-12) 

r 

and 

M  =  (D-1 3) 

9  X 

The  fact  that  x3(t)  =  ar(t)  behaves  statistically  like  a  Brownian  motion 
implies  that  the  time  derivative 

x3  =  ar  ( D-1 4 ) 

behaves  as  a  white-noise  amplitude. 


(D-l 5) 


1 

i 


The  matrix  R  =  [R.  .]  is  defined  by  the  ensemble  average  as  follows: 

~  1 J  t 

<»t(t)*j(t)>  =  f  Rjjdt 

Jo 

In  our  case  the  1  and  j  denote  the  various  pixels  of  the  detector.  No 

correlation  between  different  pixels  is  assumed,  and  it  is  also  assumed  that 


Rij 


for  each  pixel. 

Thus,  It  may  be  shown  that 


N  =  MT  R"1  M 


(D-l 6) 


(D-l 7) 


has  only  one  non-vanfshlng  element 


where 


(D-18) 


n-l 

2  w  (AlL  (D-l 9) 

63  R  J 


i 

i 


i 

) 

■wM 


l 


] 


Here  W  denotes  the  width  of  a  rectangular  target,  &  is  the  pixel 
stze,  and  AT  Is  the  contrast  at  a  target  boundary.  (Note  that  2W/6  is  the 
number  of  effective  boundary  pixels,  while  AT/6  is  the  intensity  slope  at  the 
detector  Image  boundary  due  to  finite  pixel  size) 

It  may  be  shown  that  the  quantity  R  is  related  to  noise  equivalent  temperature 

NEAT,  and  the  frame  rate  B,  as  follows: 

n  _  (NEAT)2  (D-20) 

R  -  — B - 

Defining  the  slgnal-to-nolse  ratio  SNR  by 

SNR  =  NEaT  (D-21 ) 

It  follows  that  the  matrix  element  q  may  be  written 

q  =  M  (SNR>2  B  (D-22) 


1 

* 


(! 

I! 

0 
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Steady  State  Solution  of  the  One- Dimensional  Poisson  Force  Problem 


The  steady  state  solution  of  the  Riccati  equation  is  obtained  by 
setting  the  time  derivative  of  the  variance  matrix  equal  to  zero  and 
solving  the  resulting  algebraic  equations.  MACSYMA  has  done  this  task 
for  us,  with  the  result  that  the  tracker  spatial  error  variance,  ,  is 
the  solution  of  the  following  quartic  equation: 

(Vn  q)4  +  4  D  (Vn  q)3  +  4  D2  (V^  q)2  -  S/pT  (V„  q)  -  8  D  v/p~q  =  0 

(D-23) 


while  the  other  elements  of  V  are  given  in  terms  of  V^,  as  follows: 


12 

V21 

q  vf,/2 

13  = 

V31  = 

v/p/q 

22  = 

(qV3  + 

qD  V2,  )/2  -jm 

P  s 

23 

V32 

33-  = 

(Pq  (V- 

n  q  +  2D)  Vn)/2 

(D-24) 


Consider  the  example: 
SNR  =  0.1 

B  =  15  hz 

W  =  10  ft. 

w/6  =  10 

One  f i nds 

q  *  3  ft2/sec 


( D-25) 


(D-26) 


Using  the  values  related  to  the  tank's  stochastic  dynamics  (derived 


in  a  following  section), 

P  *  0.1334  ft2/sec5 
D  *  0.025/sec 


(D-27) 
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The  quartic  equation  (D-23)  has  only  one  real  solution. 


V11 

=  0.564 

ft2 

C\J 

H 

> 

=  V21  = 

0.477 

ft2/sec 

V13 

=  V31  = 

0.211 

ft2/sec 

V22 

=  0.609 

2  2 
ft  /sec 

V23 

=  V32  = 

0.357 

ft2/s ec 

V33 

=  0.311 

2  4 

ft  /secH 

The  RMS  tracking  error  is  given  by 


(D-28) 


=  0.75  ft  (D- 29) 

i.e.,  the  RMS  tracking  error  is  75*  of  the  target's  width  for  these  para¬ 
meters  . 

Effect  of  Neglecting  D 

It  is  of  interest  to  determine  how  much  the  results  are  changed  if  the 
damping  factor  D  is  set  equal  to  zero.  One  finds 


V11 

=  2(P/q5) 

1/6 

V12 

=  V21 

2(P/q2)1/3 

V13 

=  V31 

/P7q 

V22 

-  3/p7q 

V23 

=  V32  = 

2(P2/q)1/3 

V33 

=  2(P5/q) 

1/6 

0-6 


] 

] 

(D-30)  j 

] 

] 

] 

) 

] 


i 


1 


Using  the  same  values  for  p  and  q  as  above,  one  finds 


V11 

=  0.572 

ft* 

V12 

=  V21 

=  0.491  ft2/sec 

V13 

=  V31 

=  0.211  ft2/sec: 

V22 

=  0.633  ft2/sec2  * 

V23 

=  V32 

=  0.362  ft2/sec: 

V33 

=  0.311 

ft2/sec^ 

Evidently  the  results  are  only  slightly  changed.  In  particular,  the 
change  in  the  RMS  tracking  error  estimate  is  less  than  one  percent. 

The  fact  that  Vjj  is  proportional  to  the  sixth-root  of  P,  so  that  RMS 
tracking  error  depends  upon  the  twelfth  root  of  p,  is  fortunate  since  we  probably 
can’t  estimate  P  too  accurately.  {Taking  high  roots  of  numbers  brings  them 
closer  to  1.) 
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CHOICE  OF  THE  DYNAMICAL  CONSTANTS  D  AND  P 


The  constant  D  is  obtained  from  our  model  using  data  obtained  from  a 
tank-performance  calculator.  (See  attachment).  The  calculator  gives  data 
on  the  time  to  reach  30  mph  as  a  function  of  the  tank's  horsepower  per  ton 
ratio  for  two  top  speeds--38  mph  and  63  mph. 

The  differential  equation  of  motion  may  also  be  written 

^  -  D(V  -  V) 
dt  max 

where  v  is  the  velocity  of  the  tank 
F 

and  V =  -= —  is  the  tank's  maximum  velocity  in  terms  of  the 

max  mu 

maximum  applied  force  F 

max 

Integrating  the  differential  equation,  assuming  the  tank  starts  from  rest, 
and  solving  for  D  leads  to  the  equation 

D  = 


0-8 


The  following  table  Indicates  results  obtained. 


_ Vmax  =  38  mph _ 

vtoax 

=  63  moh 

HP/TON 

t-sec 

D  — 
sec 

t-sec 

D  — 
sec 

30 

12.9 

0.121 

12.4 

0.052 

40 

9.5 

0.164 

9.0 

0.072 

50 

7.7 

0.202 

7.2 

0.090 

60 

7.7 

0.236 

6.0 

0.108 

70 

5.8 

0.268 

5.3 

0.122 

80 

5.2 

0.300 

4.7 

0.138 

The  time  t  to  obtain  30  miles  per  hour  were  obtained  by  reading  the 
tank  calculator,  while  the  values  for  D  were  obtained  from  the  above 
formula  by  setting  V=30  and  V  v  equal  to  either  38  or  63. 

iild  X 

To  obtain  P,  we  assume  that  o,  is  some  reasonable  fraction 

rr 


3(=l/2  say)  of  the  maximum  force  F( 

o 


max 
=  B  F, 


max  =  F  m  D  v  max 


therefore 


"  *  »r  <8  D  W >' 
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Applying  this  to  the  Soviet  T-62  medium  tank  (see  enclosure),  which  has  a 


580  h.p.  engine  and  weight  of  41  tons,  we  thus  get  a  HP/TON  ratio  of  14.15. 
Extrapolating  the  tabular  data  down  to  this  ratio,  and  assuming  a  top 

speed 

V  =  63  mph  =  92.4  ft/sec 
max 


leads  to 


D  =  0.025/sec 


assuming  the  tank  undergoes  a  step  acceleration  (or  deceleration)  on  the 
average  once  very  ten  seconds,  so  that 

Xf  *  0.1/sec 

one  finds  with  0  chosen  equal  to  1/2  and  V  as  above,  that 

FflaX 

P  =  0.1334  ft?/sec5 
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A  tohmn  at  T-S& 


T-62  Medium  Tank 

Anaamant:  Ona  116mm  gun.  alavafion 
ttm  16*.  dtpmtaion  mtnm  3*.  44  round*  of 
116mm  ammunition  am  carriad. 

Ona  7.62mm  machlna  gun  co-axM  with 
main  ntmamant. 

Langth:  31ft  Otn  (gun  forward) 

22ft  pin  (vuhicia) 

Haight:  7ft  10m 
Wright:  62.000  Iba. 

Manga:  260  mdaa 

6pood:  60  mg.h.  (mad) 

landing:  4ft  6M.  1 M  with  anorkal 

6ngina:  V-12  dlaaai  davaioging  660  fig. 

Craw:  4 

Width:  lift  Obi 

C/Claaranaa:  16*m 

C/Prtaaara:  11.40  pal 


V/04 await:  2ft  Tin 
Tran  ah:  OhSn 
land:  20- 140mm 


STtSFawiLai  »  dm  public  far  tfta 
*M  dma  m  Mag  1666k  and  ft  it  hthawd  w 
hwa gana Mia gmduadan m  1661/62.  Hit 


Soviet  Union 

dmalppad  from  tha  T-64  and  T-68  aariaa. 
Tha  modificationa  ineluda  a  longar  and  widar 
huA  naw  lunar  mo unrad  dlghdy  mom  to  tha 
mar.  mom  Pack  on  dm  ground,  no  and- 


tunot  hatchaa.  Tha  T-62  it  armad  with  a  naw 
118mm  gun  that  hma  hn-ataMiaad  ammu¬ 
nition  and  may  pottiMv  6a  ffttad  with  a  laaar 
rangatindar.  tha  gun  it  longar  that  that  ftfWd  to 
dm  aatfiar  T-64  and  T-86  and  hat  a  bom 
avacuator. 

Tha  T-62  can  bo  fftwd  with  auxiliary  fuai  tanba 
at  dm  roar  of  tha  vahida,  tha  T-62  can  aiao 
lay  ha  own  amoba  aeman  from  its  aahauat  pipaa 
an  aithar  aida  of  tha  hull.  Infra-wd  aguip- 
mant  la  aimdat  to  that  fttwd  to  tha  T-66.  Tha 
T-62  it  aaadv  mragnlaabla  from  tha  T-64  and 
T-66  by  ha  longar  gun  and  tha  widar  apaea 

- - am.  ---a  a - a- 


Tha  T-62  ia  aupplamandni  dm  T-64  and  T-66 
in  dm  Saviat  Army.  H  hat  alto  baan  auppiiad 
la  fgypt  and  Caachoaloyabia.  M  ia  rapaiwd 
dmt  aoma  haua  alto  baan  auppiiad  to  Poland 


AUTOMOTIVE  PERFORMANCE 


HQHtfOWlW  Ptn 


APPENDIX  E  -  DERIVATION  OF  BOBROVS KY-ZAKAI  BOUND  FOR  A  CLUTTERED  SCENE 


In  this  appendix,  the  Bobrovsky-Zakal  (B-Z)  bound  is  derived  for  the 
rectangular  target  model  Introduced  In  Section  5.1.  The  dynamic  model,  for 
both  target  and  clutter,  is 

dXT  *  Sq  dW  (target  motion)  (E-l) 

dXc  *  0  (clutter  is  at  fixed  location)  (E-2) 

where  Xj  and  Xc  are  the  target  and  clutter  location,  respectively;  q  Is  the 
targets  velocity  spectrum  amplitude,  and  W  is  a  normalized  wiener  process. 

The  measurement  equation  is 

dy  -  T(Xt,Xc)  dt  +  /r  df  (E-3) 

where  y“  Is  a  NR  •  Nj.  x  1  vector  (NR  is  the  number  of  rows  and  Nc  is  the  num¬ 
ber  of  columns)  of  pixel  measurements.  The  1th  pixel  measurement  Is  given  by 

h1(XT’XC)  *  TT  \  (XT)0-h^(Xc))  (E-4) 

+  TC  hCl  <XC> 


E-l 


where 


yv 


target  intensity, 
cl  Otter  intensity 
clutter  function  given  by 


?  /  /  P'Xe'f  rKt  (^i-)  •  dexdcr  (M) 


where 


a  *  pixel  size 
*C  ■  length  of  clutter  object 
hC  3  he*9ht  of  clutter  object 

pixel,  )  .  P  1f  **•*»  °"  »'«’  ' 
lO  otherwise 


hT^  “  target  function 


where 


•  «9 

*  7  /  / 


pixel^  (cjj.ey)  rect 


f5?) 


{^)*a 


height  of  target  object, 
length  of  target  object. 


<E-6) 


1 

I 

J 

J 

1 

1 

1 

J 

] 

*Li 

I 


In  order  to  apply  the  B-Z  bound  to  this  problem  we  first  compute 


Af  =  E  £ 


(E-7) 


where 


X  «  and  T  «  [OJ  specifies  the  dynamic  model.  Then  we  solve 

for  B  in  the  equation 

»W)''6  •  0-i)T  <CcV'  .  (§.  a)} 

♦  Ej§r  (“D1)'1  §} 


(E-8) 


C  -  specifies  the  random  process  which  drives  the  dynamic  model;  D  *  /r 

where  r  *  NEAT/Af,  Af  =  the  frame  rate,  and  D  specifies  the  noise 
level  of  the  measurement  process. 

For  the  dynamic  model  of  Interest,  equation  (E-8)  reduces  to 


Thus,  only 


=  1  E(  ahT  3h 


3X  3x 


(E-9) 


a~hT  ah 
ax  ax 


(E-10) 


need  be  evaluated  in  order  to  apply  the  B-Z  bound.  This  is  done  by  first 
noting  that 


(E— 17) 


In  order  to  evaluate  equations 
non-zero  only  for  those  pixels 


3h, 

(E-15)  through  (E-17)  we  note  that  77-  will  be 

dXr 


ftn  thp  far»pfs  Ipft 


_ 1  _  >. 


ah. 

r-r-  will  be  non-zero  only  for  those  pixels  on  the  clutter  object  left  or  right 
3xc 


edge;  that  Is, 


(E— 18) 


where 


-1  If  pixel  1  Is  on  right  edge  of  target 
1  If  pixel  1  is  on  left  edge  of  target 
0  otherwise 


and 


where 


(E-19) 


(E-20) 


1  if  pixel  1  Is  on  right  edge  of  clutter  object 
1  If  pixel  1  is  on  left  edge  of  clutter  object  (e-21) 
0  otherwise 
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For  equation  (E-15) 


I 


where 

Nc  *  number  of  pixels  on  clutter  object, 

Nn  *  number  of  pixels  on  left  and  right  clutter  edge 

Kc 

■  max  (hji  h^) 


and 


E\(Xt)  ’  "r"cl'2Nr 


nr  "c, 


EhT,  (*T>  "  Nr*, 


where 

N„  «  number  pixel  on  target  left  or  right  edge 
kT 


we  get 


(E-27) 


For  the  first  term  of  equation  (E-26)  we  observe  that  from  equation  (E-19) 

*  "V  *  P  (pixel  1  Is  on  left  or  right  edge  of  target)  (E-28) 
a 

Thus 


(F.-29) 


where  Pj  Is  the  probability  that  a  pixel  Is  on  the  left  or  right  target  edge 
and  Vj  ^Nc,Np  j  Is  as  defined  In  equation  (E-27). 

Now  for  equation  (E-16). 


(E-30) 


E-8 


where  we  have  again  used  the  assumption  that  the  target's  and  clutter's 
locations  are  independent.  Using  a  procedure  similar  to  that  used  to  derive 
equation  (E-29)  we  get 


and 


Np  *  number  of  clutter  object  pixels. 

Rather  than  evaluating  equation  (E-17),  we  observed  that  there  are  an 
equal  number  of  left  and  right  edges  on  both  the  clutter  and  target  object. 
Thus 


NR*Nr  V"c  ,  v  ,  v 

«  C  K  M  _ 

‘  J  3X^  ^3 Xy /  vXC/ 


(E-32) 


We  are  now  in  a  position  to  fill  in  the  elements  of  the  matrix  of  equa¬ 
tion'  (E-ll)  by  combining  equations  (E-12)  through  (E-14)  with  equations  (E-29) 
through  (E-32)  .  This  results  In 


E-9 


E-10 


I 


tTt*  ®t  0 

E  h  h  ■  T 

0  a, 


(E-38) 


Thus  a  solution  to  equation  (E-9)  is 


r/a^  0 


(E-39) 


Therefore  the  mean  square  error  on 


du  =  C  dw 


(E-40) 


ds*Budt  +  Vr  dZ 


given  by 


V  =  E  (iT  -TT)^  (u  -  u)  is  a  lower  bound  to  the  mean  square  error  of  the 
non-linear  problem  specified  by  equations  (E-l)  and  (E-2) . 

Thus  we  only  need  to  solve  the  matrix  Riccati  equation 


q  0  r°T  0 

-V 

0  0.  0  a. 


(E-41) 


For  the  variable  of  interest,  i.e.,  the  target  location,  we  get  for  the  mean 


square  error: 


_  (1  -  exp  (/ajq/r  t) 


(E-42) 


where  it  has  been  assumed  that  the  target  and  clutter  locations  were  known 

at  t  ■  0. 
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APPENDIX  F  PROGRAM  LISTING 


C  MAIN  ROUTINE  FOR 
C  XSTEF SIMULATING  OPTIMAT  FILTER 

C  F  (REPRESENTING  DYNAMICS),  AND  K  (REPRESENTING  OUTPUTS) 

C  ARE  TVO  SUBROUTINES  ENTERED  BY  USER  AND  LINKED  BY  USER 

CHARACTERM0  OUTFIL 
CHARACTER*5  ANS 
REAL*8  XOA , XOB 
COKI  (ON/DATA/ZT!  10000) 

DIMENSION  I LU( 7 ) , I  UP ( 7 ) 

DIMENSION  BUPPER(7), BALOWER{ 7 ) 

DIMENSION  XLOW! 7  > . XUP! 7 > , XSTEP < 7  ) , NSTEP! 7 > , XZERO! 7 ) 

com;on/specs/nfovx,nfovy,detx,dety 

COKMON/TARGET/XS I ZE , YS IZE , DELTAT , XTAR 
DIMENSION  ZMN( 7 > ,ZSIG< 7 > .ALOWER! 7 > .UPPER! 7 > 
DIMENSION  UNLWR<  7 ) ,UNUPR( 7 ) ,XX( 7  > ,STEP( 7 ) 
COKMON/VARS/U!  7  ) ,  XJA(  50000 ) ,  XJB! 50000 ) 

COMMON/NOISE /QM1 ( 7 ) , RM1 ( 10000 ) 

COHNON/ODESCR/KLEN! 7  > , SCALE ( 7  > ,OFF( 7  ) , NDIM 
DIMENSION  BOFF ( 7  ) , BSCL ( 7  ) 

DATA  BOFF ,BSCL/7*0. ,7*1 ./ 

DATA  NSTEP/7* 1 / 

WRITE (6, 1918) 

1918  FORMAT! IX , '  INPUT  SEED'  > 

READ! 5 , *  )  I  SEED 
CALL  SET_SEED( ISEED) 

C  SIMULATION  INPUTS 

WRITE! 6,135) 

135  FORMAT! IX,'  OUTPUT  FILE  NAME 7 1  ) 

READ! 5,531 )0UTF1L 
531  FORMAT! A40) 

OPEN! UNIT-20 ,NAME*OUTF IL .STATUS* ' NEW  > 

WRITE ( 6 , 881 1  ) 

8311  FORMAT! IX,*  NUMBER  OF  TRIALS’) 

READ! 5 , *  )NTRIAL 
WRITE! 6 , 2212  ) 

2212  FORMAT!  IX,’  TARGET  SIZE  X,V) 

READ(5,*)XSIZE,YSIZE 

WRITE! 6,2213) 

2213  FORMAT! IX , '  FOV#X,#Y’> 

READ! 5 , *  JNFOVX , NFOVY 
NOBS-NFOVX«NFOVY 
WRITE! 6,2214  ) 

2214  FORMAT! IX,'  PIXEL  SIZE  X,Y’) 

READ! 5 , *  )DETX ,DETY 

WRITE! 6 ,2215) 

2215  FORMAT! IX,’  DELTA-T’ > 

READ! 5,*  JDELTAT 
WRITE (6, 1000) 

1000  FORMAT! IX.'  ENTER  NUMBER  OF  STATES') 

READ! 5,*  )NDIM 
WRITE ( 6 , 1010 ) 

1010  FORMAT! IX,  'ENTER  Q  INVERSE  DIAGONAL  ELTS  ') 

READ! 5 , *  X  QM1 ( I >, 1*1 , NDIM  ) 

WRITE! 6, 1020) 

1020  FORMAT! IX,'  ENTER  R  INVERSE  ’) 

READ! 5 , *  )XRM1 

DO  2929  1*1 .NOBS 
RMl! I )-XRMl 
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2929  CONTINUE 


WRITE! 6 , 1030 ) 

1030  FORMAT! IX,'  ENTER  NUMBER  OF  POINTS,  RANGE<LOWER,UPPER> 

.  FOR  EACH  DIMENSION  *) 

READ! 5,*  XKLEN! I ) ,XLOW( I ) , XUP ! I ) , I-i , NDIM > 

WRITE! 6 , 1040  > 

1040  FORMAT! IX.’  ENTER  INITIAL  UNIFORM  DENSITV  LIMITS 
.  <LOWER , UPPER>  AND  STEP  SIZES  '  ) 

READ! 5 , * ) ( UNLWR! I > , UNUPR! I > , XSTEP! I ), I-I ,NDIM> 

WRITE! 6, 1050  > 

10S0  FORMAT! IX, *  ENTER  NUMBER  OF  ITERATIONS  ') 

READ! 5 , *  )  ITERS 
WRITE! 6 , 1060 ) 

1060  FORMAT! IX.'  WISH  TO  HAVE  STATE  PRINTOUT?  ’) 

I SWITCH-0 
READ! 5 , 1070 ) IANS 
1070  FORMAT! AI) 

IF! IANS . EQ . ' Y ’ JISWTCH-1 
WRITE! 6 , 1080 ) 

1080  FORMAT! IX, ‘  ENTER  INITIAL  STATE  VAR  VAL*’S  '  > 

READ! 5, *  X XZERO! I ) , 1-1 , NDIM) 

WRITE! 6 , 1 23  ) 

123  FORMAT! IX,1  HOW  ABOUT  A  GAUSSIAN  PRIOR?' > 

READ! 5,111  )ANS 
111  FORMAT! A5) 

C  INITIALIZE  NUMBER  OF  STEPS 

CALL  NUMSTEP! XLOW, XUP , XSTEP , NDIM, NSTEP ) 

C  START  SIMULATION  TRIALS 

DO  8033  ITRIAL-1.NTR1AL 
IF! 1TRIAL.EQ.NTRIAL  1ISWTCH-1 

C  INITIALIZE  # 

CALL  SETHALOWER, XLOW,  1, NDIM) 

CALL  SETHUPPER, XUP, 1, NDIM) 

CALL  SETUSTEP, XSTEP, 1  .NDIM) 

CALL  SETliXX, XZERO, 1, NDIM) 

C  INTIALIZE  ARRAY 

CALL  GETSCHALOWER, UPPER, SCALE, OFF) 

VAL-1 . 

C  INITIALIZE  UPPER  AND  LOWER  BOUNDARIES 

DO  50  INDX1-1 , 7 

IF! I NO XI .GT.HDIM)GOTO  40 

VAL-! UNUPR! INDX1 ) -UNLWR! INDXl ) >«VAL 

I LW( INDXl >- IF  IX! SCALE! INDXl >*( UNLWR! INDXl  l-OFF I INDXl >> )*1 
I  UP! INDXl >- 1  FIX! SCALE! INDXl  )«( UNUPR! INDXl  I-OFF! INDXl  ))  >  +  l 
GOTO  50 

40  ILW! INDXl >-l 

SCALE! INDXl >-l. 

I  UP !  INDXl  )-l 
KLEN! INDXl  )-l 
50  CONTINUE 

C  SET  U?  INITIAL  PROBABILITIES 
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VAL-l./VAL 

CALL  SETUP ( XJAi VAL  ,  ILW, IUP , ALOWER , UPPER,ANS ) 

C  START  ITERATIONS 

DO  1 00  INDX2-1, ITERS 

C  UPDATE  PROBABILITIES 

CALL  FILTER! XX ,QM1 , RM1 ) 

CALL  UPDATE! ALOWER. UPPER, STEP, ZMN.ZSIC) 

C  SWAP  ARRAYS 

CALL  SWITCH! XJA, XJB ) 

WRITE! 6,1 100 >! IQ.ZHN! IQ),ZSIC( IQ), XX! IQ > , 10-1 , ND1M ) 
WRITE! 20,1 103  1NDIM 

1103  FORMAT! IX, ‘  DIMENSION  OF  STATE  SPACE-’, 13) 

WRITE! 20, 1100)! 1Q.ZMN! IQ),ZS1S!IQ),XX!IQ),1Q-1,ND1M) 
1100  FORMAT! IX , '  VARIABLE  ',13,'  HAS  MEAN  '.F12.S, 

1  AND  SIGMA  '  ,  F 1 2 . 5 ,  '  IS  ACTALLY  ’.F12.5) 

IF! ISWTCH.EQ. 1 )CALL  PPRINTiXJA) 

C  RESET  OFFSET,  SCALE  AND  INTEGRATION  STEP  SIZE 

CALL  UP LOW! ZMN , ZS IG , BALOWER, BUPPER ) 

CALL  GETSCL! BALOWER, BUPPER, BSCL.BOFF  > 

CALL  ROOUT ! ALOWER , UP P ER , BSCL , BOF F . XOA , XOB ) 

CALL  SETHALOWER, BALOWER,  l.NDIM) 

CALL  SET11UPPER, BUPPER, 1 ,NDIM> 

CALL  SET1ISCALE, BSCL. l.NDIM) 

CALL  SETllOFF.BOFF, l.NDIM) 

CALL  SWITCH! XJA , XJB  ) 

CALL  STEPSZ!NSTEP, ALOWER, UPPER. STEP) 

100  CONTINUE 

WRITE! 6 , 1645 ) 

8888  WRITE! 20, 164S ) 

1645  FORMAT! IX,’  - - •> 

CLOSE! UNIT-20) 

STOP 

END 

C - 

C  SUBROUTINE  NUMSTEP  COMPUTES  THE  NUMBER  OF  STEPS 

SUBROUTINE  NUMSTEP! XLOW, XUP , XSTEP , ND1M , NSTEP  > 

DIMENSION  XLOW! 1 >,XUP< 1 ), XSTEP! 1  ), NSTEP! 1  ) 

DO  1  I-l.NDIM 

NSTEP! I )-INT< ! <  XUP! I )-XLOW! 1 ) )/XSTEP! I ) )♦ . 5 ) 

1  CONTINUE 

RETURN 

END 

C - - - - - - 

C  SUBROUTINE  STEPSZ  COMPUTES  THE  SIZE  OF  STEPS 

SUBROUTINE  STEPSZ! NSTEP .ALOWER , UPPER , NDIM, STEP  ) 

DIMENSION  NSTEP! 1 ),STEP! 1 >, ALOWER! 1  ), UPPER! 1 > 
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00  1  I-l.NDIM 

STEP<  I )■<  UPPER! I )-AL0VER< 1 ) >/FL0AT< NSTEP< I > ) 
1  CONTINUE 

RETURN 

END 


C  UPLOV  COMPUTES  BOUNDARY  OF  PROBABILITY  DENSITY 

SUBROUTINE  UPLOWI ZMN ,ZSIG , BLOW.BUP  ) 
DIMENSION  ZMN<7),ZSIG(7),BL0W<7>,BUP<7> 
COMMON/ ODE SCR /KLENl 7),SCALE<7),OFF<7),NDIM 

DO  1  I-l.NDIM 

BLOW< I )  — 4.*ZSIG<  I >*ZMN(  1 ) 

BUP< I >-  4 . *ZS IG< I  )*ZMN< I  ) 

1  CONTINUE 

RETURN 

END 


C  RDOUT  RESAMPLES  PROBABILITY  DENSITY 

SUBROUTINE  RDOUTI ALOWER , UPPER, BSCL, BOFF , XJA ,  XJB  ) 
DIMENSION  ALOWERI 7),UPPER(7),BSCL<7), BOFF ( 7  > , U( 7  ) 
COMMON/ JDESCR/KLEN!  7 ) , SCALE ( 7).OFF(7),NDIM 
REAL«8  X0A< 50000  l.VAL  XJINT 

REAL*8  X<JB<KLEN<  1) .  KLEN<  2  > , KLEN<  3  ) , KLEN<  4  ) , KLEN<  5  ) , 
KL EN<  6  > , KLEN<  7  )  ) 

DO  100  1 1 *1 , KLEN<  1  > 

DO  1C0  12-1 ,KLEN<2) 

DO  100  13-1 , KLEN<  3  ) 

DO  100  14-1 ,KLEN<4  ) 

DO  100  15-1 ,KLEN<  5  ) 

DO  100  1 6- 1 , KLEN<  6  1 
DO  100  17*1 ,KLEN<  7  1 
CALL  CONVRT <11,12,13,14,15.16.17, 

.  BSCL, BOFF, U> 

VAL-XO I NT<  XJA , U  ) 

XOB< 11,12,13,14,15,16,17  l-VAL 
100  CONTINUE 

RETURN 

END 


C  CONVERT  INDICES  TO  STATE  VARIABLES 

SUBROUTINE  CONVRT< II , 12,13,14,15,16,17, 
BSCL, BOFF, U) 

DIMENSION  BOFF<7),BSCL<7),U(7) 

U<  1  l-BOFF  (  1 )♦< 1 1  - 1 >/BSCL<  1  ) 

U<  2  >*BOFF<  2 )  +  < 12-1 )/BSCL<2) 

U( 3 l-BOFF <  3 )♦( 13-1 )/BSCL<3) 

U<  4  )-30FF <  4  )♦< 14-1  )/BSCL<  4) 

U( S )»BOFF<  5 )♦< 15-1 ) /BSCL <  5  1 
U<  G  J-BOFF <  6  )♦( 16-1  )/BSCL<6) 
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U( 7  )-BOFF<  7  )♦( 17-1 >/BSCL<7> 

RETURN 

END 


C  UPDATE  FILTER  AND  OUTPUTS 

SUBROUTINE  F ILTERl XX ,OMl . RMI) 
COMMON/JOESCR/KLENf 7), SCALE! 7 ),OFF< 7 ),NDIM 
DIMENSION  QM1(7>,RM1( 10000) 

DIMENSION  XX(7),ZZ< 10000), YY<  7 ) 
COMMON/DATA/ZT ( 10000) 

C  GET  NEW  STATE 

CALL  F<  XX , YY , NDIM ) 

CALL  N01V<ZZ,QM1 ,NDIM) 

CALL  AVV( YY,ZZ,XX,NDIM> 

C  GET  NEW  MEASUREMENTS 

CALL  H<  XX , ZT , NOUTS > 

CALL  N01V1ZZ.RM1 .NOUTS) 

CALL  AVV(ZT,ZZ,ZT, NOUTS) 

RETURN 

END 


C  COMPUTE  NX1  NOISE  VECTOR  A  WITH  VAR  B 

SUBROUTINE  NOIV(A.B.N) 

DIMENSION  A(N),B(N) 

DO  100  INDX-l.N 

A< INDX  )-GAUSS(0. , SQRT( 1 . /B( INDX  ) )  > 
100  CONTINUE 

RETURN 

END 


C  ADD  VECTOR  VI  TO  VZ  AND  STORE  IN  V3 

SUBROUTINE  AVV< VI , VZ , V3 , N  ) 
DIMENSION  V1(N),VZ(N),V3(N) 

DO  100  INDX-l.N 

V3( INDX )-Vl( INDX  )*VZ( INDX) 

100  CONTINUE 

RETURN 

END 


C  RESET  UPPER  AND  LOWER  BOUNDS 

SUBROUTINE  RESETMLOWER.ZMN.ZSIC.SCL  > 
DIMENSION  ALOVERl 7  )  ,ZMN< 7 ) ,ZSIG( 7  ) 
COMMON/JDESCR/KLEN<  7 ) . SCALE! 7 ) ,OFF( 7  ) , NDIM 


DO  100  INDX-l.NDIM 

ALOWER! INDX  >«ZMN( INOX >-SCL*ZSIG! 1N0X  > 
100  CONTINUE 

RETURN 

END 


- - 

C  GET  SCALE  AND  OFFSET 

SUBROUTINE  CETSCL! ALOWER, UPPER , SCALE , OFF > 

DIMENSION  ALOV/E R( 7),UPP£R(7),SCALE(7),OFF!7) 

COMMON/ JOE SCR/KLEN! 7 ) , X SCALE! 7  ) , XOFF  < 7  ) , NDIM 

DO  100  INDX-l.NDIM 

SCALE! INDX )-FLOAT( KLEN( INDX )-2 >/( UPPER! INDX l-ALOWER! INDX)) 
OFF! INDX  )>ALOWER< INDX ) 

100  CONTINUE 

RETURN 

END 

- - 

C  INITIALIZE  PRIOR  PROBABILITY  DENSITY 

SUBROUTINE  SETUP! XOA.VAL , ILW, IUP .ALOWER, UPPER, ANS ) 
DIMENSION  IUP! 7 ) , ILW! 7 ) 

CHARACTER*5  ANS 

DIMENSION  ALOWER! 7),UPPERI7) 

COt'.MON/JDESCR/KLEN! 7 ) .SCALE! 7 ) .OFF! 7 ) .NDIM 
REAL*8  XJAIKLEN! 1  ) , KLEN! 2 > , KLEN! 3  ) , KLEN! 4 ) . KLEN! 5 > . 

.  KLEN! 6  )  ,KLEN! 7  >  ) 

DO  100  1 1 “1 , KLEN! 1  ) 

DO  100  1 2* 1 , KLEN! 2  > 

DO  100  1 3“ 1 , KLEN! 3 ) 

DO  100  I4-1.KLENU) 

DO  100  15*1 .KLEN! 5  ) 

DO  100  16* I , KLEN! 6 ) 

DO  100  17*1  .KLEN! 7  > 

XOA! II, 12, 13, 14, 15, 16, 17) *0. 

100  CONTINUE 

IFIANS.EQ. 1  YES  1 .OR . ANS . EQ . * Y  ') GOTO  2001 

C  IF  UNIFORM 

DO  200  1 1 «ILW! 1 ), IUP! 1 > 

DO  200  I2-ILW! 2  ) , IUP! 2 ) 

DO  200  I3*ILW( 3  ) , IUP! 3  > 

DO  200  1 4*ILW<  4  ) , IUP ( 4 ) 

DO  200  I5-ILV! 5 ) , IUP! 5 ) 

DO  200  I6-1LW! 6  ) , IUP! 6 ) 

DO  200  1 7*ILW( 7 ) , IUP!  7 ) 

XJA! 11,12,13,14,15,16,17 )*VAL 
200  CONTINUE 

GOTO  1000 

C  IF  GAUSSIAN 

2001  DO  201  11*1 *KLENi 1 ) 
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201 


DO  201  I2-1.KLENI2) 

DO  201  13-1 ,KLEN<3> 

OO  201  14-1 , KLEN!  4 ) 

00  201  I5-1,KL£N<5) 

DO  201  16-1 , KLEN( 6 ) 

DO  201  17-1.KLEN17) 

XDA!  11(12,13,14,15,16,17  J-GDENSITY! 11, 12, 13. 14, IS. 16.17. I LW. I UP, 

ALOWER, UPPER  > 

CONTINUE 


1000  CONTINUE 

RETURN 

END 


C  COMPUTE  GAUSSIAN  PROBABILITY  DENSITY 

REAL *8  FUNCTION  GDENSITYI 11, 12, 13, 14, 15. 16. 17. I LW.t UP .ALOWER, UPPER) 
DIMENSION  ALOWER! 7),UPPER(7) 

COMMON/UDESCR/KLEN! 7 ) .  SCALE!  7  )  .OFF!  7  > , NDIM 
DIMENSION  INDEX!7),ILW(7),IUP(7) 

REAL'S  0 , ZDEN 

INDEX! 1  )-It 
INDEX! 21-12 
INDEX! 31-13 
INDEX! 41-14 
INDEX! 51-15 
INDEX! 6  >  —  1 6 
INDEX! 71-17 

D-l . 

DO  10  I-1.ND1M 

XMU-FLOAT!  ILW!  1 1  I'FLOAT!  1UP<  l  )-HV<  1 1)12. 

SIGMA-FLOAT! IUP! I l-ILWI I) 1/8. 

X-FLOAT! INDEX! I ) ) 

D-D*ZDEN! X.XMU, SIGMA)*! ! KLEN! I )-2 )/! UPPER! I l-ALOWERI I ) ) > 

10  CONTINUE 


GDENSITY-D 

RETURN 

END 


C . . - - - - - 

C  COMPUTE  GAUSSIAN  PROBABILITY 

REAL'S  FUNCTION  ZDEN! X , XMU, S IGMA > 

REAL'S  D 

D-! X-XMU ) /SIGMA 

ZDEN-DEXP l-.5*D*D>7! SIGMA 'SORT! 2 . *3 . 1 4 ) > 

RETURN 

END 


C - - - 

C  PRINT  RESULTS 

SUBROUTINE  PPRINTIXDA) 

COMMON 70 DESCR/KLEN! 7  ) , SCALE! 7  >  .OFF! 7  ) , NDIM 
REAL'S  X0A1KLEN! 1 ) .KLEN! 2  ) , KLEN! 3 > .KLEN! 4  ) , 
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KLEN< 5 > ,KLEN! 6 > ,KLEN< 7 )  )  • 

WRITE! 6, 1000)! <(<((( XOA( 11.12, 13, 14, IS, 16,17), 11*1, KLEN(l)), 

12-1, KLEN<  2  )  > , 

1 3* 1 , KL EN( 3  > ) , 

14*1 , KLEN! 4  > ) , 

15*1 , KLEN( 5 ) ) , 

16*1 , KLEN( 6  >  > , 

17-1 , KLEN( 7  )  > 


WRITE ( 20, 1000 )(ii!l!!X0A!Il,I2,I3,14,I5,I6,I7),Il-l,KLENi  1  >), 

12*1 , KLEN! 2)1, 

13- 1 , KLEN! 3 ) )  • 

1 4- 1 ,KLEN! 4  > ) , 

15- 1 , KLEN! 5 ) ) , 

16- 1 ,KIEN<  6 ) > , 

17- 1 ,  KLEN! 7  >  ) 


1000  FORMAT! IX, 5D 12. 5) 
RETURN 
END 

C - 


C  SWAP  ARRAYS 


SUBROUTINE  SWITCH! XOA.XOB  ) 

COMMON/ JDESCR/KLENi 7 ) , SCALE! 7  > ,OFF<  7 ) , NDIM 
REAL*8  TEMP 

REAL*8  XJA! KLEN! 1 ) , KLEN! 2 } , KLEN! 3  ) , KLENl 4  ) , 
.  KLEN!5),KLEN(6),KLEN(7)> 

REAL*8  XOBiKLENi 1 ) , KLEN! 2  ) , KLEN! 3  ) , KLEN! 4  ) , 
.  KLEN! 5  > , KLEN! 6  ) , KLEN! 7  )  ) 


DO  100  I1-1,KLEN!  1  > 

DO  100  12*1 , KLEN! 2  ? 

DO  100  1 3- 1 , KLEN! 3  ) 

DO  100  14-1 , KLEN! 4  > 

DO  100  15-1 , KLEN! 5  ) 

DO  100  16-1 , KLEN! 6  ) 

DO  100  17- 1 , KLEN! 7  1 
TEMP-XOA! II ,12,13,14,15,16,17) 

XJA! 1 1, 12, 13, 1 4, 15, 16, 1 7 l-XJB!  1 1,12, 13, 14, 15, 16, 17) 
XOB! II , 12, 13, 14. 15,16, I7)-TEMP 
100  CONTINUE 


RETURN 

END 


C  COMPUTE  IMAGE 


SUBROUTINE  Ht V , HOUT , NOUT  ) 

DIMENSION  Vi  7), HOUT! 10000) 
COMMON/SPECS/NFOVX , NFOVY . DETX ,DETY 
COMMON/ TARGET/ XS I 2E , YS IZE , DELTAT, XTAR 
NOUT* NFOVX* NFOVY 
ICNT-0 


C  TARGET  LOCATION 


XTAR-V! 1 ) 

YTAR— YSIZE/2.+NFOVY*DETY/2. 
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DO  100  0-l.NFOVY 
DO  100  I-l.NFOVX 
ICNT-1CNT*1 

C  CONVERT  PIXELS  1,0  TO  IMAGE  UNITS 

YO-FLOAT! J-l >*DETY-FLOAT( NFOVY-1 >*DETY/2. 
XI-FLOAT( 1-1 l-DETX -FLOAT! NFOVX-1 l-DETX/2. 

C  PIXEL  BLUR  FACTOR 

XFACT-DETOUT! XTAR ,X1,DETX,XS1ZE ) 

C  PIXEL  OUTPUT 

VAL«XFACT*DELTAT 
HOUT< ICNTl-VAL 
100  CONTINUE 

RETURN 

END 

C - 

C  COMPUTE  PIXEL  BLUR  FACTOR 

FUNCTION  DETOUT! X.PICLOC, BLUR. SIZE) 

C  TARGET  LOCATION  IS  ZERO 
XX-PICLOC-X 

C  CHECK  TO  SEE  IF  TARGET  COVERS  ANY  OF  THIS  PIXEL 
ABSXX-ABS!XX> 

IF ! ABSXX .CT .  .6*<  SIZE+BLUR ) >COTO  50 

C  IS  PIXEL  ON  LEFT  EOGE  OF  TARGET? 

XLEFT-XX*! SIZE-BLUR >/2. 

IF! XLEFT.GT.0. 1GOTO  60 

C  PIXEL  IS  ON  LEFT  EDGE  OF  TARGET  SO  US  TRI.  FUNCTION 
C  TO  COMPUTE  PIXEL  OUTPUT 

ARG"! XL EFT) /BLUR 
DETOUT-TRIIARG) 

RETURN 

C  PIXEL  NOT  ON  LEFT  EDGE  SO  CHECK  RIGHT  EDGE 

60  XRIGHT*XX-! SIZE-BLUR )/2. 

IF ! XR1GHT . LT.0. )GOTO  70 

C  PIXEL  IS  ON  RIGHT  EDGE  OF  TARGET  SO  US  TRI.  FUNCTION 
C  TO  COMPUTE  PIXEL  OUTPUT 

ARGa<  XRIGHT  l/BLUR 
OETOUT«TR1!ARG) 

RETURN 

C  PIXEL  IN  MIDDLE  REGION  OF  TARGET 
70  DETOUT- 1. 
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RETURN 


C  PIXEL  NOT  ON  TARGET 

S0  DE TOUT-0. 
RETURN 
END 


C - 

C  TRIANGLE  FUNCTION 

FUNCTION  TRI<  X ) 
V-ABS(X) 

IF(V.CT.l. )GOTO  1 
TRI-l.-V 
RETURN 
1  TRI-0. 

RETURN 

END 


C - 

C  STATE  TRANSITION  FUNCTION 

SUBROUTINE  F( V.FOUT.NSTATE ) 

DIMENSION  V< 1  ),FOUT< 1  ) 

NSTATE-1 
FOUTI 1  )-V< I > 

RETURN 

END 

C - - - 

C  MW 

C  THIS  ROUTINE  WILL  COMPUTE  THE  SCALAR  TRANPOSEIA )*B*C 
C  WHERE  B  IS  DIAGONAL.  A,C  ARE  VECTORS  OF  DIMENSION  N 

SUBROUTINE  MVV< A, B ,C , N . OUT > 

DIMENSION  A<N).B<N).C(N> 

REAL-8  OUT 
OUT-0. 

C  8  CONTAINS  ONLV  THE  DIAGONAL  ELEMENTS 
DO  100  INDX-1 ,N 

OUT-OUT+AI INDX >*B( INDX  >«C( INDX  > 

100  CONTINUE 

RETURN 

END 

C - - - - - 

C  DVV 

C  THIS  ROUTINE  WILL  COMPUTE  A-B  WHERE 
C  A.B  ARE  N  DIMENSIONAL  VECTORS 

SUBROUTINE  DVV(A,B,C,N) 

DIMENSION  A( N  )  ,B<  N ) ,C( N ) 

DO  100  INDX-1, N 

C< INDX  )-A< INDX  )  —  B £ INDX) 

100  CONTINUE 
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RETURN 

END 


C - 

C  XJINT 

C  THIS  ROUTINE  WILL  INTERPOLATE  THE  XO 
C  PROBABILITY  DENSITY  FUNCTION  TO  GET 
C  TRUE  VALUE.  INPUT  IS  A  VECTOR  U.  OUTPUT 
C  IS  THE  SCALAR  VALUE. 

REAL*8  FUNCTION  XJINT(XJ.U) 

REAL*0  XO , A, B , XOGRAB 

DIMENSION  XJI 50000 >,ISUBSC(7),DX(7),U(7> 
COMMON/JDESCR/KLEN! 7  > , SCALE! 7 > ,OFF! 7 > , NDIM 

C  GET  INDICES 

DO  100  INDX-1,7 

C  HANDLE  HIGHER  DIMENSION  NOT  USED 
IC-1 

I F ( INDX.CT.NDIM )GOTO  50 

C  CTS  OF  U  FOR  INDEX 

PLACE-SCALE! INDX  )•< U< INDX  >-OFF( INDX  )  > 

IC-IF IX( PLACE )  +  l 

C  ERROR  CHECKING 

IF( ( IC-l.LT. 1  ).OR.< IC+l.GT.KLEN! INDX))>GOTO  500 

C  GET  DELTA  FOR  INTERPOLATION 

DX< INDX  )-( PLACE-FLOAT! IC-1 > ) 

50  ISUBSC! INDXI-IC 

100  CONTINUE 

C  NOW  DO  INTERPOLATION 

A-XJGRAB! XJ , ISUBSC ) 

DO  200  INDX1-1 , NDIM 

C  INCREMENT  LOOKUP  INDEX 

ISUBSC! IN0X1  )■ I SUBSC! INDX1  >  +  l 
B-XJCRAB(XJ, ISUBSC) 

C  INTERPOLATE  BY  ADDING  IN  DERIVATIVE 

A-A+! B-A)*0X( INOXI  ) 

C  RESTORE  SUBSCRIPT 

ISUBSC! INOXI >-lSUBSC! INDXI >-l 
2 00  CONTINUE 

C  OUTPUT  RESULT 

IF! A. LT.0. )A-0. 
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XJINT-A 

RETURN 


C  ERROR  CHECK 

5 00  CONTINUE 

XJINT=0. 

RETURN 

END 


C - - - 

C  XJGRAB 

C  FUNCTION  TO  GET  AN  ELEMENT  FROM  THE 
C  XJ  ARRAY,  ACCESSED  BY  THE  COEFFICIENTS 
C  IN  THE  ISUSSC  ARRAY 

REAl*8  FUNCTION  XJGRAB! XJ , I SUBSC > 

DIMENSION  ISUBSC1 7 ) 

COMMON/JDESCR/KL  EN! 7  > , SCALE! 7  )  ,OFF< 7  )  ,NDIM 
REAL*8  XJIKLENU  > ,KLEN< 2 > ,KLEN< 3 > , KLEN! 4 > , 
KLEN<S),KLEN(6),KLEN(7>) 

XJGRAB-XJ( ISUBSCl 1 ) . I SUBSC! 2 > , I SUBSC! 3 > . ISUBSCt 4 ) , 
ISUBSCl 5  ) , ISUBSC! 6 > ,  ISUBSC! 7 > ) 

RETURN 

END 


C  XOPUT 

C  SUBROUTINE  TO  PUT  AN  ELEMENT  IN  XJ  ARRAY 

SUBROUTINE  XJPUT< XJ . VAL , ISUBSC J 

COMMON/ JDESCR/KLEN!  7  >, SCALE!  7  )  ,OFF(  7  ) , NDIM 

DIMENSION  ISUBSCl 7) 

REAL*8  VAL 

REAL*8  XJ<  KLEN( 1 ) ,KLEN< 2 ) ,KLEN( 3 )  ,KLEN( 4  ) , 

.  KLEN(5),KLEN(G),KLEN(7>) 

XJ( ISUBSCl I ) , ISUBSCl Z ) , I SUBSC! 3  )  .ISUBSCl 4  )  .ISUBSCl 5  ) , 
.  ISUBSCl 6), ISUBSCl 7 )>-VAL 

RETURN 
END 


C - 

C  UPDATE 

C  UPDATE  PASSES  THROUGH  AN  ITERATION  OF 
C  THE  FILTER  EQUATION 

C  IS  SETS  UP  XJB  FROM  XJA  AND  PRODUCES  MEANS  AND  SIGMAS 

SUBROUTINE  UPDATE! ALOWER , UPPER . STEP 
.  ,ZMN,ZS1C> 

REAL *8  XJA , XJB , VAL ,ZNUM, SUMI NT, FUNUM, FX1 ,FX2 
COMMON/JDESCR/KLENI 7 ) , SCALE! 7 ) , OFF <  7  > , NDIM 
COMMON/ VARS/U! 7  )  ,XJA< 60000 ) ,XJB! 50000  > 

COMMON / I VARA/ I VAR 
DIMENSION  JSUBSC! 7  ) , IONE! 7  ) 

DIMENSION  ALOWER! 7),UPPER!7 > .STEP! 7 ) , ZMN! 7 ) 
DIMENSION  ZS IG! 7  > 

DATA  IONE/7*l / 

EXTERNAL  FUNUM, FX1 ,FX2 

C  INITIALIZE  JSUBSC, AND  U 
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CALL  SET1  ( 0SU8SC , IONE ,1,7) 

CALL  SETl! U .ALOWER , 1,7) 

100  LDIM-1 

C  CALCULATE  NUMERATOR  OF  FILTER  EXPRESSION 

ZNUM«SUM1NTIFUNUN,ND1M, ALOWER. UPPER, STEP) 
CALL  X0PUT( XJB ,ZNUM, JSUBSC ) 

C  INCREMENT  SUBSCRIPT  AND  U 

2 00  0SUBSC<  LDIM)-JSUBSC( LDIM  HI 

U< LDIM  )-U< LDIM  HI. /SCALE! LDIM) 

C  FINISH  THIS  DIMENSION? 

IF IOSUBSC! LDIM ).LE.KLEN! LDIM >)GOTO  IN 

C  YES-INCREMENT  DIMENSION-RESET  U 

JSUBSCf LDIM  )■ 1 
U( LDIM  ) ■ALOWER! LDIM) 

LDIM-LDIMM 

C  ARE  YOU  FINISHED  WITH  ALL  DIMENSIONS? 

IFiLDIM.LE .  NDIM)SOTO  200 

C  NORMALIZE  PROBABILITY  DENSITY 

CALL  NORMLZE < ALOWER , UPPER , STEP , XOB ) 

C  YES-GET  MEAN  AND  SIGMA 
C  TAKE  MEANS  AND  SIGMAS  OF  STATE  VARIABLES 

00  300  INOX- l, NO IM 

C  SET  I VAR  TO  THE  BVARIABLE  INDEX  TO  BE  OPERATED  ON 
I VAR- INOX 

C  FX1  CONSISTS  OF  U*0!U> 

ZVR-SUMINT<FX2,NDIM. ALOWER, UPPER, STEP) 

ZMN< INDX )-SUMINT(FXl . NDIM, ALOWER , 

.  UPPER, STEP) 


C  CALCULATE  SIGMA 

I F C ZVR-ZMN!  INDX )*ZMN! INDXI.LT.0. > WRITE! 6 , 121Z > 
1212  FORMAT! IX,1  NEG  SIGMA  SHIT  ’) 

ZSIG! INDX )*SQRT! ( ABS! ZVR-ZMN! INDX )*ZMN! INDX  >  >  >  > 
300  CONTINUE 

RETURN 

END 

C - 

C  NORMALIZE  PROBABILITY  DENSITY 

SUBROUTINE  NORMLZE!  ALOWER, UPPER, STEP, XM> 
COMMON/JDESCR/KLEN! 7 ) , SCALE! 7 )  .OFF! 7  ) , NDIM 
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REAL *8  XJBI KLEN< 1 > .KLENl 2 ) .KLENl 3  )  .KLENl 4 > , 
KLEN(5),KLEN(6),KLEN(7>) 

DIMENSION  UPPER! 7  )  ,ALOWER( 7),STEPI7) 

REAL-8  DENSITV ,SUM,HTERM 
EXTERNAL  DENSITV 

DO  2 000  II-l.KLEN(l) 

DO  2000  12-1 .KLENl 2  ) 

DO  2000  13-1. KLENl 3) 

DO  2000  14-1. KLENl 4) 

DO  2000  15-1. KLENl 5) 

00  2000  16-1, KLENl 6) 

DO  2000  17-1. KLENl 7) 

XJBI 1 1,12, 13. 14. IS. 16, 17  l-XOBI 11,12.13.14,15.16,17) 

.  -HTERMI 1 1. 12. 13. 14, IS, 16. 17, ND1M, SCALE. OFF) 

2000  CONTINUE 

SUM-SUMINTI DENS ITV.NDIM.ALOWER, UPPER. STEP) 

DO  200  11-1 .KLENl 1  ) 

DO  2 00  12-1 .KLENl 2  ) 

DO  200  13-1 .KLENl 3 ) 

DO  200  14-1, KLENl 4) 

DO  200  15-1, KLENl 6) 

DO  2 00  16-1 .KLENl 6  ) 

00  200  17-1, KLENl 7) 

XJBI 11,12,13,14,15,16,17  )-XJBI 11,12,13,14,15,16,17 )/SUM 
200  CONTINUE 

RETURN 

END 


C  CURRENT  CONDITIONAL  PROBABILITY  DENSITV 

REAL-8  FUNCTION  DENSITY! X) 

DIMENSION  XI 7) 

REAL-8  XJA.XJB.XJINT 

COMMON /VARS /U l 7 ) .XJAI 50000 ) , XJBI 50000) 
DENSITY-XJINTIXJB.X) 

RETURN 

END 


PERFORM  NUMERICAL  INTEGRATION  IN  UP  TO  10  DIMENSIONS 
CF  ARBITRARY  FUNCTION  . 

INPUTS! 

FUNC  FUNCTION  TO  BE  INTEGRATED.  MUST  BE  DEFINED  AS 
AN  EXTERNAL  IN  USERS  CALLING  ROUTINE. 

NDIM  DIMENSIONS  OF  SPACE  TO  BE  INTERGATED  OVER. 

FROM  NDIM  OIMENIONAL  VECTOR  WHICH  SPECIFIES  THE  LOWER 
INTERCATION  LIMIT  FOR  EACH  DIMENSION. 

TO  NDIM  DIMENSIONAL  VECTOR  WHICH  SPECIFIES  THE  UPPER 

INTEGRATION  LIMIT  FOR  EACH  DIMENSION. 

STEP  NDIM  DIMENISIONAL  VECTOR  WHICH  SPECIFICES 

THE  INTERGATION  STEP  SIZE  FOR  EACH  DIMENSION. 

REAL-8  FUNCTION  SUMINTIFUNC, NDIM, FROM, TO, STEP ) 

REAL-8  FUNC, TOTAL, ZERO 
DIMENSION  FROM! 10), TO! 10), STEP! 10) 

DIMENSION  XI 10),ZERO(0*10> 

DIMENSION  TOTAL 10! 10) 
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DATA  ZERO/11*#./ 

C  INITIAL  TOTAL  TO  ZERO 

CALL  S£T( TOTAL , ZERO ,0, NDIM ) 

C  INITIALIZE  X-  VARIBLE  OF  INTEGRATION. 

CALL  SET1 ( X , FROM, 1 iNDIM) 
l 00  I  DIM- 1 

C  COMPUTE  INTERGRAN 

TOTAL ! 0 )«FUNC( X  > 

C  INCREMENT  X 

200  X<  IDIMI-X!  IDIMI+STEP!  ID1M) 

TOTAL! IOIM)-TOTAL< IDIMI+TOTAL! IDIM-1 >*STEP< IDIMI 

C  RESET  LOWER  DIMENSION  SUM 

TOTAL! IDIM-1 >«0. 

IF !  X!  IDIMI.LT.  TO!  IDIMDCOTO  100 

C  RE INITIL IZE  LOWER  DIMENSION 

X! IDIMI-FROM! IDIM) 

C  GO  ONE  DIMENSION  HIGHER 

IDIM-IDIM+1 

C  FINSHED? 

IF! IDIM.LE.NDIMIGOTO  200 

C  GET  SUM 

SUM1NT-TOTAL1NDIM) 

RETURN 

END 

C - 

C  SET  VECTOR  X  TO  VECTOR  V,  DOUBLE  PRECISION 

SUBROUTINE  SET! X , Y .NSTART, NSTOP > 

REAL*8  X! NSTART i NSTOP  ) , V! NSTART i NSTOP  ) 

DO  100  I -NSTART, NSTOP 
X!  I  >-Y!  I  ) 

100  CONTINUE 

RETURN 

END 

C - 

C  SET  VECTOR  X  TO  VECTOR  Y,  SINGLE  PRECISION 
SUBROUTINE  SETliX, Y, NSTART, NSTOP > 
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DIMENSION  X<  NSTART :NSTOP ) , VI NSTART i NSTOP ) 


DO  1 00  1 -NSTART, NSTOP 
X< I >•¥<  I  > 

100  CONTINUE 

RETURN 

END 


C  FUNUM- INTEGRAND  FUNCTION  FOR  FILTER'S  NUMERATOR 

REAL*8  FUNCTION  FUNUM(V) 

REAL*8  XJA,X0B,T4.T8.DD.DTl,DT2,XJINT,T2 
DIMENSION  V( 7 ) , HOUT ( 10000 ) . FOUT( 7  >, TEMPI  10000) 
COMMON/VARS/UI 7 ) , X0A<  60000 ) , XJB< 50000 > 

COMMON/ DATA /ZT( 10000) 

COMMON/NOI SE/QM1 < 7 ) , RM1 < 10000  > 

C  COMPUTE  SECOND  TERM 

CALL  F( V.FOUT.NSTATE  ) 

CALL  DVV<  U ,FOUT , UMV , NSTATE  ) 

CALL  MVV< UMV, QM1 , UMV, NSTATE, T2 ) 

DT2-T2 

DD-DEXP<-.S«(DT2)) 

C  LOOK  UP  OLD  PROBABILITY 

T4«XJINT(XJA,V) 

T5-T4*DD 

FUNUIt-TS 

RETURN 

END 


REAL*8  FUNCTION  HTERM< 1 1 , 1 2, 13, 14 ,  IS , 16 , 17 , NDIM, SCALE .OFF  ) 
REAL*8  DTI 

COMMON/DATA/ZT ( 10000) 

COMMON/NOISE /GM1 ( 7 ) , RMl ( 10000  > 

DIMENSION  SCALE <  7),OFF(7),U(7) 

DIMENSION  TEMPI  10000), HOUTI 10000) 

CALL  CONVRT < 11,12,13,14,15,16,17, SCALE .OFF , U ) 

CALL  H<  U , HOUT , NOUTS  ) 

CALL  DVVIZT, HOUT, TEMP, NOUTS) 

CALL  MVVITEMP, RMl, TEMP, NOUTS, TI  ) 

DT1-T1 

HTERM-DEXPI - . S*0T1  ) 

RETURN 

END 


C 


REAL*8  FUNCTION  FXKV) 

DIMENSION  V<7> 

REAL*8  XOA , XJB , X J INT 

COMMON /VARS / U ( 7 ) , X J A ( 60000 ) , X 0 B ( 60000 ) 

COMMON/ I VARA/ 1 VAR 

FXl-XJINTI XOB ,V  )*V( I VAR ) 
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RETURN 

END 


REAL *8  FUNCTION  FXZIV) 

DIMENSION  V(7) 

REAL *8  XOA.XOB.XJINT 

coi;mon/vars/u( 7  > , xoai 500001 ,xob( 50000) 

COMMON/ I VARA/ I VAR 

FX2-XJINT! XJB , V )  *V! IVAR)*V( IVAR) 

RETURN 

END 


THIS  FUNCTION  USE  THE  REJECTION-ACCEPTANCE 
METHOD  TO  GENERATE  A  GAUSSIAN! XMEAN, SIGMA > 
RANDOM  VARIBLE. 

FUNCTION  GAUSS! XMEAN, SIGMA) 

COMMON /SEED/ISEED 

GENERATE  TWO  UNIFORMS, 1) 

1  Ul-RAN! ISEED) 

U2-RAN! ISEED) 

X— ALOGIU1  ) 

IFIU2.GT.EXPI-IX-1. >«*2/2. > )GOTO  J 


C  X  HAS  BEEN  ACCEPTED.  GENERATE 
C  RANDOM  SIGN  AND  ATTACH  TO  X. 

P-RAN! I  ) 

SIGN-1. 

IF!P.GT..5)SIGN— 1. 

X-SICN*X 

C 

C  X  IS  NOW  GAUSS! 0,1).  CONVERT  TO 
C  GAUSS! XMEAN, SIGMA). 

GAUSS-SICMA*X*XMEAN 

RETURN 

END 

C - 

SUBROUTINE  SET-SEED! I ) 

CHARACTER'S  TM 
COMMON/SEED/ 1  SEED 
DATA  I  SEED/99999/ 

C  IF  I<-  0  SET  SEED  FROM  CLOCK 

C  OTHERWISE  SEED  IS  SET  TO  USUER  SUPPLIED  VALUE  I 
IF!I.CT.0)GOTO  10 
CALL  TIME(TM) 

I l-ICHAR! TM! 9 1 S ) > 

I SEED-2*! 11*11*1 00 ) ♦ I 
RETURN 
C 

10  ISEED-I 
RETURN 
END 
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